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2 . 0 ABSTRACT 


The  scientific  computer  program  SIGNIR  was  written  to  provide 
predictions  of  the  infrared  emission  from  axisymmetric  turbojet, 
turbofan  and  turboshaft  engine  exhaust  system  hot  surfaces.  This 
manual  provides  the  program  user  with  all  information  necessary 
to  obtain  predictions  of  the  infrared  signature  of  any  selected 
exhaust  system.  The  program  conducts  its  own  internal  flow  analysis 
and  thermal  balance  to  provide  the  system  surface  temperature  dis- 
tribution. Both  surface  radiation,  including  multiple  surface  re- 
flections, and  special  surface  cooling  methods  are  employed  by  the 
program.  The  radiant  intensity  of  the  emitted  energy,  attenuated 
by  the  atmosphere,  is  provided  for  selected  ranges  and  system  off- 
axis  angles  for  infrared  wavelengths  of  from  one  to  15  microns  at 
0.05  micron  increments.  Also,  selected  infrared  radiant  intensity 
bandwidth  ourmations  may  be  requested  which  can  include  the  detector 
response  characteristics. 

The  ph/sical  exhaust  system  selected  is  translated  into  a 
mathematical  program  model  by  the  user.  All  required  internal 
program  decisions  are  keyed  from  the  model  input  data.  This 
manual  provides  the  user  with  a description  of  the  program,  the 
analytical  techniques  employed,  the  methods  required  to  model  the 
selected  exhaust  system,  program  input  data  requirements  and  output 
data  available.  The  user  will  find  that  little  Knowledge  is  re- 
quired, of  the  program's  internal  operations  to  be  capable  of  utiliz- 
ing program  SIGNIR  for  exhaust  system  signature  predictions. 
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3.0  PROGRAM  DESCRIPTION 


The  description  of  the  infrared  signature  prediction  program  * 

(SIGNER)  is  provided  by  the  following  sections  which  cover  a general 
description  of  the  program,  the  program  computational  procedure  and 
the  analytical  methods  utilized  in  obtaining  the  solution.  A sec- 
tion covering  special  data  computation  is  provided  to  aid  the  program 
user. 


3.1  GENERAL 

The  program  SIGNIR  is  a digital  computer  program  written  for 
the  purpose  of  providing  predictions  of  the  hot  metal  infrared  emis- 
sion from  aircraft  engine  exhaust  systems.  The  program  is  designed 
to  be  applicable  to  axi- symmetric  turbojet,  turbofan,  or  turboshaft 
engine  exhaust  systems.  It  predicts  the  spectral  intensity  of  the 
radiant  energy  emitted  from  exhaust  system  hot  metal  in  the  wave- 
length band  of  from  one  to  15  microns.  In  general,  the  information 
required  by  the  program  is  as  follows: 

. Exhaust  system  physical  characteristics. 

. Engine  operating  conditions. 

. Special  surface  cooling  flow  conditions. 

. Exhaust  system  surface  properties. 


The  predictions  provided  by  the  program  for  the  combination  of  a 
selected  maximum  of  30  engine  off-axis  angles,  five  detector  ranges 
and  three  detector  altitudes  are: 

. Spectral  radiant  intensity 

. Selected  radiant  intensity  bandwidth  summation 


Optional  exhaust  system  information  which  can  be  requested  from  the 
program  is: 

. Internal  fluid  flow  properties. 

. Sprface  boundary  layer  data. 

. Internal  and  external  geometric  radiation  view  factors. 

. Surface  temperature  distribution. 

. Thrust  information. 
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Program  SIGNIR  is  made  up  of  71  subprograms.  This  program  is  written 
in  Fortran  IV  for  the  IBM  360  system.  A program  flow  diagram,  a program 
listing  and  a basic  task  description  for  each  of  the  subroutines  are 
provided  in  the  Appendices . All  necessary  computer  control  information 
required  to  operate  this  program  is  discussed  in  section  7. 

3.2  CONPUTATICNAL  PROCEDURE 

The  process  of  computing  the  intensity  of  the  infrared  energy  emitted 
from  aircraft  exhaust  system  hot  metal  requires  a knowledge  of  the  metallic 
surfaces  temperature,  emissive  properties  and  their  relative  physical 
positioning  within  the  system.  From  this  information,  the  radiant  energy 
emitted  and  reflected  from  the  exhaust  system  to  a point  in  space  can  be 
determined.  The  computational  procedure  established  by  this  program  for  the 
prediction  of  exhaust  system  hot  metal  infrared  radiation  is  shown  by  major 
subjects  in  the  general  block  diagram  presented  on  figure  1. 

The  input  data  required  by  the  program  provides  a physical  model  (nodal 
system)  of  the  exhaust  system  configuration,  the  properties  and  conditions  of 
the  exhaust  gases  entering  the  system.  The  program  undertakes  a one -dimensional, 
isentropic,  compressible  flw  analysis  which  generates  the  axial  distribution 
of  exhaust  gas  flow  conditions  and  properties  within  the  system.  Where  two 
gas  streams  at  different  temperatures  flow  within  a single  passage,  compound 
flow  considerations  are  utilized.  The  program  adjusts  the  flow  at  the  nozzle 
throat  to  balance  with  the  ambient  pressure  or  a chocked  condition  dependent  on 
the  system's  entrance  flow  data. 

With  the  internal  flew  properties  defined,  a surface  boundary  layer  analysis 
is  conducted  to  obtain  the  coefficients  of  convection  and  friction  along  each 
surface  within  the  system.  This  analysis  is  primarily  dependent  on  the 
surface's  axial  static  pressure  distribution  and  the  upstream  boundary  layer 
conditions.  An  area  weighted  average  coefficient  of  convection  for  each  of  the 
model's  surface  nodes  (isothermal  regions  established  by  the  input  data)  is 
computed  for  the  system  thermal  analysis.  The  axial  force  on  each  surface, 
internal  to  the  system,  is  computed  utilizing  both  the  frictional  and  pressure- 
area  forces. 


^■Reduced  by  ASD 
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The  next  step  taken  by  the  program  i~  the  calculation  of  the 
temperatures  along  any  surfaces  cooled  by  ; ran spi  r at  i on , film  or 
convection-film  cooling  methods.  The  program  accepts  a fixed 
coolant  flow  rate  or  computes  a system  coolant  flow  balance  when 
coolant  source  conditions  are  defined.  The  calculation  for  these 
cooling  methods  are  dependent  on  the  physical  design  of  the  wall 
itself  and  the  conditions  of  the  gas  flow  adjacent  to  it.  The 
axial  surface  temperature  distribution  is  an  area  weighted  average 
across  each  of  the  nodes  involved  along  the  cooled  surface. 

The  calculation  of  internal  system  geometric  radiation  view 
factors  is  the  next  step  in  the  programs  computational  procedure. 

One  of  the  more  difficult  factors  in  radiation  calculations  is 
determining  the  fraction  of  energy  leaving  one  given  surface  that 
will  be  intercepted  by  another  given  surface,  or  a portion  of  that 
surface.  Geometric  view  factors  provide  this  description.  This 
program  computes  view  factors  for  each  node  taking  advantage  of 
the  symmetry  of  the  system  and  considering  each  node  as  a frustum 
of  a right  circular  cone,  a ../linder  or  a disk.  View  fa  tors  are 
computed  for  each  node's  view  of  all  other  nodes  within  the  sys- 
tem including  itself.  The  computation  procedure  accounts  for  sha- 
dowing when  the  view  between  two  nodes  is  partially  obscured  by 
some  Intervening  surface. 

With  the  geometric  radiation  View  factors  computed  and  system 
surface  emissivities  available  from  input  data,  the  overall  radia- 
tion interchange  factor  (script  F)  is  computed.  This  factor,  used 
in  the  radiation  calculations,  accounts  for  the  view  between  two 
surfaces,  the  non-black  nature  of  these  surfaces  and  the  fact  the 
multiple  reflections  between  all  surfaces  are  occurring.  Two  sur- 
faces, which  geometrically  can  not  see  each  other,  can  interchange 
energy  through  energy  reflection  from  other  surfaces  within  the 
system.  The  computational  procedure  considers  all  surfaces  to  be 
diffused  surfaces,  i.e.  the  emitted  and  reflected  energy  departing 
a surface  has  the  same  distribution  as  that  of  a black  body  which 
follows  Lambert's  law  of  cosines. 

At  this  point  in  the  computational  procedure,  all  essential 
information  necessary  to  conduct  a system  internal  heat  balance 
has  been  acquired.  The  steady  state  thermal  analysis  considers 
energy  radiation,  convection  and  conduction  The  conduction  be- 
tween nodes  along  the  same  surface  is  neglected  since  the  thermal 
gradient  along  the  walls  are  expected  to  be  low.  Radiation  from 
the  hot  gas  has  also  been  neglected  since  this  energy  contribution 
to  a surface  is  small  when  compared  to  the  convection  from  the  hot 
' gas  and  or  the  radiation  from  exhaust  system  surfaces  at  near  the 
same  temperatures.  The  programs  thermal  analysis  provides  the  tem- 
perature of  each  of  the  surface  nodes  within  the  exhaust  system. 
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of  Jcco“plished  by  the  Program  is  the  determination 

2 radiation  interchange  factors  between  a detector  fixed  in 
P the  nodes  established  for  the  exhaust  system.  The  view 

factors  are  computed  between  the  detector  and  those  nodes  which  can 

factor^*  by  the  ?eJeC^°r‘  ^ tUrn’  overan  radiation  interchange 
fetors  are  computed  which  consider  the  multiple  energy  reflection 

which  occur  between  nodes,  the  emissive  properties  of  the  surfaces 

-^vnn?ew°me^riC  View  factors  as  was  done  Previously  with  the  in- 
ternal interchange  factors. 

+he  overall  interchange  factors  and  system  node  tempera- 
> r^ant  energy  passing  from  the  system  hot  metal  to  the 
6 °rr  utilizing  the  Stefan-Boltzmann  law.  In  turn 

the  energy  distribution  is  computed  considering  the  radiant  energy  ’ 

^°^e+^LStTlbuI;ed  in  accordance  with  Planck's  law.  This  energy  dis- 

intpn^+n  1h-P^aCi?  ^ m°re  ca!mon  terms  of  spectral  radiSt 
ntensity  which  eliminates  any  consideration  of  the  area  of  energy 
collection  by  the  detector.  energy 
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3.3  ANALYTICAL  METHODS  INCORPORATED 

The  analytical  method  utilized  by  program  SIGNIR  Is  Pre®«n^®^ 
in  this  section  under  the  general  topics  of  exhaust  syi stem  gas  flow, 
surface  cooling,  radiation,  internal  heat  balance  and  radiation  emis- 
sion. By  presenting  the  analytical  methods  in  this  ^uence,  the 
progression  of  information  is  similar  to  that  presented  by  the  com- 
putational procedure  of  section  3*2. 

3.3.1  Exhaust  System  System  Gas  Flow 

An  internal  gas  flow  analysis  is  conducted  by  SIGNIR  under  the 
control  of  RIMAIN.  The  primary  effort  is  to  determine  the  distribu- 
Son  of  gas^t  propertie  s Ld  the  surface  boundary  layer  convec- 
tion coefficients.  The  following  sections  cover  the  analytical  pro- 
cedures utiH  zed  in  computing  the  compressible  flow  and  surface 
boundary  layer. 

3. 3.1.1  Compressible  Fluid  Flow 

The  sub-program  COFLOW  establishes  the  flow  properties  within 
a compound  stream  based  on  a one-dimensional,  isentropic  !U^" 

lysis  The  method  employed  is  based  on  the  analysis  presented  in 
inference  1 which  assumes  that  the  primary  and  leeoaA^r  8tr*“s^ 
a mixed  turbofan  exhaust  system  can  be  handled  separateiy.  I ry 

has  been  found  to  agree  well  with  experimental  f 

vines  an  array  or  compound-flow  properties  associated  witn  now  area 
^r  a gi venT comb i nation  of  flow  rates.  The  fluid  static  pressure  is 
chosen  to  be  constant  across  adjacent  streams.  If  a single  stre 
configuration  is  employed,  the  same  analytical  procedure  is  employed 
utilizing  only  the  single  stream. 

The  nomenclature  used  within  this  section  is  as  follows: 

„ 2 

A - flow  area,  ft. 

CF  - flow  coefficient  2 

a - acceleration  of  gravity,  ft. /sec.  . 

j - an  integer  (corresponding  to  the  subscript  j; 

k - ratio  of  specific  heats 
& - mass  flow  rate,  lb^/sec. 

M - flow  Mach  number 

n . - number  of  compound  streams 
p - . static  pressure,  lb ft. 2 
PI  total  pressure,  lbf/ft.2  2 

AP  - increment  in  static  pressure,  lb^/ft. 
r - gas  constant,  lb^ft./lb^^ 

T - static  temperature, q°R 

Tm  - total  temperature,  R 

V - flow  velocity,  ft. /sec. 

v - compound-flow  indicator 

f - density,  lb^/ft.3 
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subscripts: 

amh  - ambient  conditions 

ex  - nozzle  exit 

i - an  integer,  stream  number 

i - an  integer,  relates  physical  conditions 

and  flow  properties  at  the  same  static 
pressure  conditions 

The  sub -program  COFLOW  computes  an  array  of  fluid  properties  com- 
patible with  values  of  static  pressure  decremented  from  the  total  p 
sure  provided  to  the  ambient  static  pressure,  i.e. 


P * PT  - J-AP 

j 


(1) 


The  stream  properties  for  a perfect  gas  are  computed  by  the 
following  isentropic  flow  relationships: 


• j <£nlt*  if-”*  - 1 . 


(2) 


Tij  - TriJ  [/  +XXj)  ■ 


(3) 


Vij  --  Mij  & ri,j 


W 


/°£  ; = Pj  / UJ  . 

LJ  J 


(5) 


Aitj  • nii  / ( ft', j Aij)  ■ 


(6) 
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A compound  flow  indicator,  is  used  to  determine  if 
choking  occurs  at  a static  pressure  value  which  is  greater 
than  ambient  pressure  (choked  convergent  nozzle).  Compound- 
flow  is  subsonic-  when  /&  > 0,  sonic  at  /$  - 0 and  supersonic 
when  < 0.  The  indicator  is  evaluated  by  the  equation 
given  in  reference  1 as 


& ■ f -/j 

i s / * w 


This  shows  that  at  the  compound- choke  condition,  some  streams 
can  be  supersonic  while  others  are  subsonic.  For  a single 
stream,  n = 1,  = 0 when  the  stream  mach  number  is  one.  The 

sub-program  COFLCW  is  limited  to  two  streams  (i  = 2).  Once  the 
array  of  stream  properties  are  formed,  the  total  flow  area, 


I*. 


at  the  nozzle  exit  is  determined  from  the  program  at 


/d  * O , //"  P > 


or  at 


P*  P* 


svnb  / 


'/*  /3  ro 


The  total  flow  area  of  the  array  is  adjusted  by  a flow  coefficient 


/ A 


which  is  an  adjustment  toward  two-dimensional  flow  considerations. 
3. 3.1.2  Surface  Bounds 


The  conditions  within  the  compressible  turbulent  boundary 
layer  which  exists  along  each  of  the  exhaust  system  surfaces  is 
computed  by  sub-program  TURBLT.  This  is  a segment  of  an  existing 


- ^''fyravTirif,  ^ 
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computer  program  written  by  A.  H.  Ybarra  of  Vought  Aeronautics 
Company  of  LTV  Aerospace  Corporation.  The  anaS^ti cal  methods 
employed  are  basically  those  presented  by  reference  2 in  which 

cralb3  ^ iS  expressed  in  tenns  of  the  momentum  inte- 

grai  and  moment- of -momentum  integral,  simplified  by  a Mager-type 

ZTi Ybarra  Lo^ratefL 
empirical  modi1  ication  to  the  moment -of -momentum  equation  which 
improved  the  methods  agreement  with  test  data. 

1 list  of  ^ 


e 

FI 

F2 

H 


h 

h 

M 


c 


m 


r 

S 


T 

U 

u 

v 


J 


■x 

Y 


<f 1 


- transformation  parameter 

local  skin  friction  coefficient,  dimensionless 

- transformation  parameter 

specific  heat  at  constant  pressure,  Btu/lb  -°R 

- exact  differential  m 

- exponential  base,  e = 2. 718291328,  dimensionless 

- empirical  factor 

- empirical  factor 

- boundary  layer  shape  factor,  dimensionless 
(incompressible  transformed  denoted  bv  sub- 
script tr) 

- enthalpy,  Btu/lb 

convection  heat  transfer  coefficient,  Btu/hr.ft2oR 

- Mach  number,  dimensionless 

- power  law  profile  exponent,  dimensionless 

- pressure,  lb  /in? 

- Prandtl  number,  dimensionless 

- Radius  of  symmetry,  ft. 

Reynolds  number  based  on  momentum  thickness, 
dimensionless 

- gas  constant,  lb  ft. /lb  °R 

- recovery  factor,  dimensionless 

- enthalpy  ratio  factor,  S = (h  Jh  )-l.O,  dimen- 
-  temperature,  °R 

- transformed  velocity  along  surface,  ft/sec. 

- physical  velocity  along  surface,  ft/sec. 
transformed  velocity  normal  to  surface,  ft/sec. 

- physical  velocity  normal  to  surface,  ft/sec. 

- transformed  coordinate  along  surface,  ft. 

- physical  coordinate  along  surface,  ft. 

- transformed  coordinate  normal  to  surface  ft. 

- physical  coordinate  normal  to  surface  ft. 
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A 

A* 

S 

6* 

G 

T > 

P 

v 

r 

a 

Subscripts 

e 

1 
o 

ref 

s 

T 

tr 

w 

T 

2 


ratio  of  fv y hiW.p,  dunrupiou’ioss 

- transformed  boundary  layer  thickness,  ft. 

- transferred  boundary  layer  displacement 
thickness,  ft. 

- physical  boundary  layer  velocity  thickness,  ft. 

- physical  boundary  layer  displacement  thickness, 
ft. 

- boundary  layer  momentum  thickness,  ft. 

- absolute  viscosity,  lb/sec-ft. 

- kinematic  viscosity,  fT^/sec. 

- density,  lb  /ft3 

- shear  stress,  lb^/ft-sec^ 

- transformed  velocity  potential,  l/sec. 

- velocity  potential,  l/sec. 

- partial  derivative  syinbol 


local  invsicid  flow,  at  edge  of  boundary  layer 
incompressible 

free stream  stagnation  conditions 
at  Eckert  reference  enthalpy  conditions 
local  stagnation  conditions 

total  conditions,  unsubscripted  indicates  static 

conditions 

transformed 

at  the  surface,  wall 

~4-  nw/j  ^ '"•enent 

previous  value,  at  start  of  current  increment 


For  steady,  compressible, turbulent  flow  the  boundary  layer 
equations,  in  which  the  flow  variable  appear  as  time-averaged  quanti- 
ties, are  expressed  by  the  continuity  and  momentum  equations  for  axi- 
symmetric  system  as 

d(j°u)  JfrwJ  d?  . o (12) 

T* ~ d*  E d* 


and 


4 /°LT 


du 

9Y 


dP  ay 
dx  &Y 


(13) 
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These  two  partial  differential  equations,  in  their  present  form, 
,,,be  solved  simultaneously.  To  simplify  and  reduce  the 

ty  inteeration)  a Mager-type  transformation  is  performed 
which  transforms  the  physical  coordinates  as 


'J 


b(-x)  dx 


and 


(14) 


V ‘ ccx)  f £ d? 

4 * 


where 


T±L 


t>M  - (Tt.  / 


CM  - (T*/Tn) 


>/l 


(15) 


(16) 

(17) 


The  velocities  can  be  replaced  through  the  definition  of  a 
stream  function,  as 


& 

ad  _ _ /o  tr 

7% 


(18) 


and  the  transformed  velocities  by  the  relations 


u 


ad 

9Y 


and 


. _ ad 
ax 


(19) 


V 


LI  — J 


w 
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Applying  equations  14,  15  and  18  to  the  boundary  layer 
equations  12  and  13  result  in  the  following  transformed  con- 
tinuity and  momentum  equation: 


€¥  + <2¥  * o (2°) 

dX  dY 


9X  dY 


o >r* 


(21) 


where  the  enthalpy  term  is  defined  as 

S s (hs/he)  - t.o  (22) 

where  h is  the  local  stagnation  enthalpy, 
s 

The  boundary  conditions  applicable  to  the  equations  20  and  21 

are: 


ZA^OJ'O 

VCX,o)‘0 


no  slip  a l wail. 

no  mass  added  or  removed  from  boundary 
layer  at  wall. 


S(Xo)*$uOC)  enthalpy  ratio  at  wall  depends  only  on 
local  recovery  temperature. 


(23) 


// inii  S » o 


local  stagnation  temperature  approaches  free- 
stream  total  temperature. 


Unit  II , j /y) boundary  layer  velocity  approaches  external 
€ velocity. 


Some  correlations  between  transformed  and  physical  parameters 
which  result  from  the  Mager-type  transformation  are  as  follows: 

Longitudinal  velocity,  a m /Tro  f*  jj  « [ccx)]  IS . 

Velocity  potential,  ^ * 5^. 


(24) 

(25) 
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Boundary  Layer  Shape  Factor,  /y  « ^ ^ J 

r±L- 

Momentum  thickness,  <£?  <git_  (7^-0  /7i)  2<'rw^  (27) 

Sr-  / 

Displacement  thickness,  y (77t>/T<r)  zCr-7)  (28) 

“ &tr(7ro/7Z:  )2(Y’'1) . 


The  combining  and  manipulation  of  equations  20  and  21  produce 
the  equation 


J -'[zsftt-wj*  ‘§p.[v(u.-ir)]  + Zgd&.v) 

+ „ ili  r - - tJ  €l¥ 

* u’  Sic  =?r‘ 

(29) 

• 

The  transformed  momentum  thickness  &tr, 
placement  thickness,  Z*are  defined -as 

and  the  transfromed  dis- 

Otr  • f (H  )('  ~ H ) dyr 

f « *-»€/  » *-«  / 

(30) 

**  ■/('  'K)dY 

(3D 

and  the  transformed  boundary  layer  shape  factor  follows  as 


Mir  - A*/Otr 


(32) 


Integrating  equation  29  with  respect  to  Y between  the  limits  Y=0 
and  Y=  A (where  A is  a distance  normal  to  the  surface  sufficiently 
large  such  that  the  conditions  S=0  and  v = are  both  satisfied), 
applying  the  chain  rule  and  equations  30  and  31 > the  following  equa- 
tion form  can  be  obtained  : 


d&tr 

dX 


1 

Li  dX  L 


Z&tr  > 


(33) 


i.  ^ ..las--,  <ir^i 
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The  fundamental  definition  of  the  boundary  layer  shear  stress 
is  expressed  by 


By  utlizing  equations  l4,  24,  3 2 and  34  and  a reference  temperature 
technique  such  that 


t (35) 


equation  33  becomes 

d@tr  + &tr  dJJtfp 
dX  U*  dXL 


(36) 


which  is  the  form  of  the  transformed  momentum  integral  equation 
presented  in  reference  2. 

The  moment-of-momentum  integral  equation  is  derived  in  a similar 
manner.  The  differential  form  the  momentum  equation  29  is  multiplied 
by  the  parameter  Y and  the  results  inturn  integrated  with  respect  to 
Y from  Y=0  to  A . The  resulting  integral  equation  obtained  is. 


dJx  ~ [“tr  (Nth  * (Ntr  " l)BtrfS  d Y 


B(jjtr-l) 


Utr 


■r+l)&tr  J 


SYdY 


+ (l - /)  / jjjp  1 / Tnrf  ) 

@tr~  ( / ( 7fo  / ££* 

- (tJtr  -l)(tJtr+l)  I f2z*£  I f Y ,/ Y J 

&tr  ( 7Z/1  TroJ  J %/  ( "Z ) 


(37) 
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The  Ludwig-Tillman  skin-friction  relationships  for  incompressible 
turbulent  flow,  both  convenient  and  valid  over  a wide  range  of  pressure 
gradients,  is  given  as 


c/ 


V „ -G.2£o6 


(38) 


The  Reynolds  number  based  on  momentum  thickness  is  defined  as 

<S>/V  - (39) 

With  the  application  of  the  reference  enthalpy  concept,  equations  2b 
and  27,  and  the  relations 


(40) 

(41) 


equation  38  can  be  put  in  tne  lorm 


4l 

z 


s 


7i/ 


O,  /Z3  (7k/ Tne f ) ( Mnrf/M*) 

(Ue&r  /U) 


m 


This  is  the  expression  of  skin  coefficient  for  compressible  flow  in  terms 
of  the  transformed  velocity,  momentum  thickness  and  reference  tem- 
perature (and  viscosity). 


The  integrals  which  appear  in  the  transformed  momentum  and 
moment-of -momentum  equations,  36  and  37,  shall  be  reduced  to  a 
usable  form.  The  recovery  factor  is  defined  as 


* « 
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for  a turbulent  boundary  layer.  For  a Prandtl  number  of  unity, 
the  recovery  temperature  will  be  equal  to  the  freestream  stagna- 
tion temperature.  Therefore,  for  adiabatic  flow,  the  value  of  S 
in  equation  22  would  bt  zero  and  the  integrals  will  be  zero.  To 
evaluate  the  integrals  at  an  arbitrary  surface  temperature,  the 
Crocco  relation  is  employed  together  with  a power  law  assumption 
for  the  velocity  profiles,  i.e. 


/is  *~  /i*> 

ho  — hu> 


U_ 


U_ 

ZXr 


(WO 


and 


(45) 


Since  the  purpose  of  the  Mager-type  transformation  is  to 
render  the  boundary  layer  equations  analogus  to  the  incompres- 
sible form,  the  transformed  and  incompressible  shape  factor  are 
near  equal.  Using  this  information,  the  integrals  are  evaluated 
as 


/J 

J dY 


Sdi&tr  ' l)Ui8tr 


and 


/ 


A 

SY  dY 


m 


(1*7) 


Substituting  equations  22,  24,  38  and  46  into  equation  36,  the 
momentum  equation  becomes 


d&tr  , f T&fjf  Q I - &&  fz  + to?  Ur]  (48) 

dX  [rnJ(sJ  ct  Sx Ls  * h.  . 


the  same  procedure,  the  term  of  equation  37  containing  the 
;egration  in  Y reduces  as 


2 

(t-Ji  + 1) 


Y/Vt' 


Ul  Z + 4 Ui  - / 
(Fi+l)(Fi+3)  * 


(49) 


e integral  contained  in  the  last  term  of  equation  37  is  not  so 
sily  simplified.  An  empirical  relation , 


£.003075  -O.OQ33SZ 
Cf/z 


(50) 


was  established  from  the  data  presented  in  reference  2 in  terms  of 
the  variables  and  Cf . Once  the  analytical  model  was  complete, 
further  empirical  corrections  in  the  last  term  of  equation  37 
found  necessary  to  correlate  with  measured  test  data.  These  corre- 
lation factors  take  the  form  of 

Ft  * &/0  + ir  (51) 


(52) 


FZ  * /.O 


uihrn  ^ - O 
dX 


= /.O  + (4*4?)'  cuh<m  < O 
l dX  / J dX 


\ 


I 

'A 

I 

] 


with  the  further  restriction  that 


(F/)f F2)  * /.O 


(53) 
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By  employing  those  methods  required  in  establishing  equation  48 
and  applying  the  empirical,  parameters  presented,  equation  37 
reduces  to  the  form 


d±k  a _ ti:fa+/)Y4-/)[lJh«,^j\ i _/_  dif* 
dX  2 l u*  cfX 


+ A QJi+O* (T,*f\  Cf 

&tr  VTtoJ  2 


- (Pl)(Fi) <yLdMd(7h t)  CcoscnsUi-oosss.z)  (:*> 


In  subprogram  TURBLT,  the  two  transformed  boundary  layer 
equations  48  and  54  are  solved  simultaneously  using  an  iterative 
finite  difference  numerical  integration  procedure,  specifically, 
the  process  employed  divides  the  surface  being  analyzed  into  »nw*n 
increments.  For  each  surface  increment,  a value  for  H.  and  Q±r 
is  assumed  at  the  downstream  end  of  that  increment  andthen,  using 
the  average  value  for  all  pertinent  parameters  over  the  surface 
increment,  the  two  equations  48  and  54  are  solved.  The  results 
are  compared  with  the  initial  estimates  and  an  iteration  carried 
out  to  converge  on  a solution. 

The  program  is  provided  with  the  properties  of  the  gas  stream 
flowing  adjacent  to  the  surface  and  the  initial  (upstream)  values 
of  momentum  thickness ; 0 i and  incompressible  shape  factor  H.  (these 
two  values  must  be  input  by  the  program  user  which  requires  some  es- 
timate of  the  boundary  layer  conditions  entering  the  exhaust  system). 
With  a surface  increment’s  upstream  (2)  boundary  layer  condition,  the 
program  makes  an  assumption  as  to  the  downstream  (1)  growth.  The  con- 
ditions and  properties  values  c-er  the  increment  are  averaged,  equa- 
tions 42,  51  and  53  are  evaluated  and,  inturn,  equations  48  and  54  are 
solved  simultaneously  utilizing  the  equations 


Qtr,  ' (X -Xt)  (55) 


on 
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and 


^ - %g(x;-x*) . 


(56) 


Th®. ^ tial  estimate  is  corrected,  the  solution  process  continued 
until  the  parameters  converge.  In turn,  the  parameters  for  the  fol- 
lowing increments  are  solved  in  a systematic  order. 

From  the  solution  obtained,  local  physical  boundary  layer  para- 
meters along  the  surface  are  computed  by  the  program.  These  are 

boundary  layer  thickness, 

s * * D ] 


local  skin  friction  coefficient. 


Cf  ” 6. 


-I.5L  V-i 


<3.-2.  C.8 


local  convection  heat  transfer  »•<««*• 


(58) 


- 

■v^ - . 

..  _ ~/3 - .Z(s,zsr  -Z2s)  M* ] (59) 


3.3.2  Exhaust  System  Surface  Cooling 

Three  method  for  cooling  exhaust  surfaces  have  been  incor- 
porated into  the  program  SIGNIR.  They  are  film,  convection- film 
and  transpiration  cooling.  These  cooling  methods  were  selected 
as  the  more  efficient  means  of  cooling  surfaces.  They  give  the 
program  latitude  in  the  types  of  systems  which  this  program  is 
applicable.  The  following  section  presents  the  analytical  pro- 
cedures employed  for  each  of  the  surface  cooling  methods. 
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3. 3. 2. I Film  Cooling 


The  gaseous  film  cooling  method  incorporated  by  subprogram 
FILMCL  is  an  empirical  correlation  of  various  investigator's  data 
presented  by  reference  3.  The  cooling  configuration  utilized  by 
the  program  is  specifically  for  tangential  injection  of  the  coolant. 
The  method  proposed  by  reference  4 was  coupled  with  this  film  cool- 
ing method  to  be  capable  of  handling  multiple  slot  configurations. 
The  n<*HHiclaturr  uisd  within  thi*  Mirtion  is  as  follows : 


i 


A 

Cd 

CP 

f(v) 

g 

h 

i 

K 

m 

M 

• 

m 

n 

P 


R 

s 

St 

T 

T 

r 

T 

T 

UA 

V 

x 


r 

A 


2 

area,  ft. 

flow  discharge,  nondimensional 

specific  heat  at  constant  pressure,  Btu/lb.  . °R 

m 

velocity  mismatch  factor,  nondimensional 

p 

gravitational  acceleration,  ft. /hr. 

convective  heat  transfer  coefficient,  Btu/hr.ft.2  °R 

an  integer,  nondimensional 

p 

pressure  loss  coefficient,  lbf/ft. 

total  number  of  cooling  slots,  nondimensional 

heat  capacity  ratio,  (oVCp)c/(/3 VCp)g 

mass  flow  rate,  lb  /hr. 

m 

flow  exponent,  nondimensional 

p 

static  pressure,  lb  /ft. 

® 2 

total  pressure,  lb  /ft. 

m 

gas  constant,  lb.. ft. /lb  . °R 
i m 

slot  height,  ft. 

Stanton  number  (h  //JVC  ),  nondimensional 

^ ^ o 

temperature  or  static  temperature,  R 

recovery  temperature,  °R 

total  temperature,  °R 

overall  rate  of  heat  transfer,  Btu/hr.  °R 
velocity,  ft. /hr. 

distance  downstream  from  slot,  ft. 

ratio  of  specific  heats 

the  change  in  a property,  nondimensional 
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^ - film  effectiveness 

^std  ” density  at  standard  temperature  and  pressure,  lb^/fl 
0-  - ratio  of  density  to  density  at  standard  temperature 

and  pressure,  nondimensional 
Subscripts: 

c - coolant 

f - fluid  film 

f’  - adjacent  to  film 

g - mainstream  gas 

s - slot 

t - total 

w - wall 

1 - coolant  source 

2 - heat  source 

The  empirical  equation  for  film  effectiveness  proposed  by 
reference  3 is 


-eiereiice  o Ty 

y-  sm]} 


where 


--  (Tr3  - T„)/(Tr,  - Tre  ) , 

./f'W  - / + Orf  src-ian  - /J  > 

(Vg  ) /-5’  “ O 

* l V6j  9 

uihan  ^ > / 

Vc  * 

u>A<rn  —!  £ / 

\/c 

These  equations  were  converted  into  a form  where  the  independent 
variables  are  those'  available  to  the  program.  Specific  terms  of 
equation  60  were  altered  in  form  fs  shown. 


hg  % Cd  As 
Cpc  rhc 


3 


(60) 

(61) 

(62) 


S6 

M 


(63) 
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7mm 


rl  -aXes  “he  .'orm 

v .-  •>.  ‘ fV  (ft) 

where  T_,  is  the  gas  film  temperature  resulting  from  the  upstream 
slot  flow.  The  term 


Ji  , _±£_ . £4.  JL 

rtc\P9J 


(66) 


from  equation  62,  is  converted  as  shown  which  assumes  isentropic 
flow  from  a common  coolant  plenum  to  the  slot  whose  discharge  static 
pressure  adjusts  to  equal  the  mainstream  gas  static  pressure.  Com- 
bining the  forgoing  equations,  the  equation  for  wall  temperature  be- 
comes 


where  the  term  f(v)  is  defined  by  equations  62  and  66. 

The  computational  process  conducted  by  FIIMCL  starts  with  the 
upstream  slot,  utilizing  the  local  gas  stream  recovery  temperature 
for  Tfl,  compute  the  wall  temperature  from  the  slot  to  the  end  of 
the  surface  being  cooled.  This  is  accomplished  in  incremental  steps 
using  the  local  gas  stream  properties  along  the  surface.  Progressing 
to  the  next  slot,  the  gas  film  temperature  (Tfl)  along  the  remaining 
surface  is  assigned  the  values  of  the  wall  temperature  (T  ) computed 
from  the  previous  slot.  In  this  manner,  upstream  slot  coolant  films 
affect  all  surface  temperatures  downstream.  The  procedure  is  con- 
tinued progressively  along  the  cooled  surface  producing  the  wall  tem- 
perature distribution  as  demonstrated  in  figure  2.  Inturn,  an  area- 
weighted-average  temperature  is  computed  for  each  of  the  nodes  assigned 
along  this  cooled  surface. 
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Mainstream  Gas  Flow  at  Temperature  T 


g 


T 

g 


t 

Temperature 


T 

c 


Surface  Temperature 

Temperature  of  Fluid  (adjacent 

to  coolant  film) 


(b)  Multiple  slot  film  cooled  surface  temperaturf  distribution 


FIGURE  2:  Multiple  Slot  Film  Cooling 
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A coolant  flow  balance  is  incorporated  in  this  subprogram 
which  provides  the  coolant  flow  properties  required  by  the  eoua- 
tions  fob  and  (>7.  The  program  user  is  provided  the  option  of ‘de- 
fining either  the  coolant  flow  rate  or  its  source.  If  flow  rate 
is  defined,  coolant  supply  temperature  oust  also  be  provided. 

The  slots  are  assumed  to  be  supplied  from  a common  plenum  such 
that  the  driving  total  pressure  to  each  slot  is  identical! 

The  coolant  flow  from  each  coolant  slot  is  a function  of  both  the 
slot  flow  area  and  the  gas  stream  local  static  pressure.  The 
coolant  total  pressure  and  individual  slot  coolant  flow  rate  is 
computed  using  an  iteration  process.  A form  of  the  continuity 
equation,  using  isentropic  considerations,  is  shown  in  terms  of 
the  dependent  variables,  slot  flow  rate  and  total  pressure. 


The  selected  total  coolant  flow  rate  is 

rr\ 

-■  y Me;  . (69) 

L */ 


The  program  iterates  about  total  pressure  and  flow  rate  to  obtain 
a solution  which  provides  the  individual  slot  coolant  flow  rates. 

If  the  option  is  selected  to  compute  the  coolant  flow  rate 
based  on  coolant  source  temperature  and  pressure,  an  extension  of 
the  forgoing  process  is  required.  This  computational  procedure  in- 
cludes the  supply  system  pressure  loss  and  heat  transfer.  The  cool- 
ant supply  system  pressure  loss  is  computed  by  the  equation  form 


0 ’ APr  - A?  fm)n 


(70) 


where  the  system  characteristics  K and  n are  required  input  (see 
discussion  in  section  3«^).  The  program  solves  this  equation  in 
the  form 

(pT'-  K'ZfTr,  +Trc)]/n  (?1) 
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The  coolant  pieman  temperature  is  evaluated  by  the  equation 


% - rTZ  -C-T^-m)  e(^/‘r‘rnc) 


(72) 


where  temperature  of  the  heat  source  (T_J  and  coolant  supply 
overall  heat  transfer  coefficient  (UA)  is  required  input  data 
(see  discussion  in  section  3.4).  The  combination  of  equations 
68,  69,  71  and  72  are  solved  by  an  iteration  process  which  pro- 
vides the  resulting  slot  flow  rates  and  coolant  plenum  tempera- 
ture and  pressure  required  in  determining  the  wall  temperature. 

3. 3*2. 2 Convection-Film  Cooling 

The  convection- film  cooling  method  computational  procedure 
incorporated  by  subprogram  CONFIM,  developed  by  VAC,  utilizes 
the  empirical  film  cooling  correlation  presented  by  reference  3, 
the  multiple  slot  method  suggested  by  reference  4,  and  the  char- 
acteristics of  plate-fin  heat  exchanger  materials  as  presented 
in  references  5 and  6.  • 

The  nomenclature  used  within  this  section  is  as  follows: 

„ 2 

A - area,  ft. 


d 

d 

e 

exp 

fF 

f(v) 

g 

G 


h 

K 

L 


- specific  heat  at  constant  pressure,  Btu/lb  . °R 

m 

- the  differential  of  a variable 

- hydraulic  diameter,  ft. 

- exponential  (e) 

- Fanning  friction  factor,  nondimensional 

- velocity  mismatch  factor,  nondimensional 

O 

- acceleration  of  gravity,  ft. /hr. 

- gas  mass  flow  velocity  through  heat  exchanger  core, 
lb  /hr. ft.2 

Cl 

- convective  heat  transfer  coefficient,  Btu/hr.ft.2  °R 

p 

- pressure  loss  coefficient,  lbf/ft. 

- length  of  convection-film  cooled  shingle,  ft. 
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- mass  flow  rate,  lb ^/hr. 
static  pressure,  lb  /ft.1- 

- °randtl  number,  nondimensional 

2 

- total  pressure,  lb  /ft. 

Cl 

- heat  flow  rate,  Btu/hr. 

- gas  constant,  Ib,ft./lb.  °R 

I m 

- Reynolds  number,  nondimensional 

- slot  height,  ft. 

- temperature  or  static  temperature,  °R 

- total  temperature,  °R 

- overall  rate  of  heat  transfer,  Btu/hr.  °R 

- distance  downstream  of  slot,  ft. 

- density,  Ib^/ft.^ 

- density  at  standard  temperature  and  pressure,  lb ^/ft. 


SUBSCRIPTS: 


- coolant 

- coolant  discharge 

- fluid  film 

- adjacent  to  film 

- mainstream  gas 

- heat  exchanger 


x - wall  location 

x+ax  - incremented  wall  location 

Figure  3 (a)  and  (b)  show  sketches  of  sectioned  convection/film 
cooling  panel  configurations  analysis.  Referring  to  the  element 
shown  on  figure  3 (c),  the  heat  transfer  from  the  fluid  film  to  the 
wall  is  equaled  to  that  transferred  from  the  wall  to  the  coolant  when 
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Mainstream  Gas  Flow 


Coolant  Supply 
(From  Plenum) 


(b)  Parallel  Flow  Configuration 


x x + dx 


(c)  Element  of  Cooling  Panel,  Energy  Balance  for 
Counterflow  Scheme 

FIGURE  3:  Convection- Film  Cooling  Scheme 
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conduction  along  the  wall  is  ignored,  i.e, 

- h,Au  - (hA)hK(T„  - -Ttc)  . 

The  wall  temperature  can  be  solved  as 

- m)[-rif/ChA)nt+  7*/(VwJ 

where 

^i/4^  * £ >/ (hA)fo  + I /(hfAui)J  # 


(73) 


(74) 


(75) 


The  heat  flow  from  the  heat  exchanger  to  plenum  coolant 
is  considered  low;  small  difference  in  temperature  with  a re- 
latively low  plenum  side  coefficient  of  convection.  The  heat 
gained  by  the  heat  exchanger  coolant,  due  to  the  predominance 
of  the  heat  flow  from  the  film,  is  found  to  be 


* 


• . — /hr  / I / 1 / . 

WcCp  a /ic  = ( J ( Uf  - J CJK 


(76) 


When  substituting  equation  74  into  this  equation  and  rearranging, 
the  equation  becomes 


OJA) 

nhcCpL 


G/X 


(77) 
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Holding  the  film  temperature  constant  over  the  incremental  dis- 
tance AX  end  integrating,  the  heat  exchanger  coolant  temperature 
is  given  by 


^ci,‘Ax=  <2pL)J  (78) 


vrtiere  the  average  film  temperature  over  the  increment  is  computed 
by 


T*f  * l 


(79) 


The  film  temperature  is  evaluated  by  methods  presented  in 
section  3. 3*2.1  where,  for  this  case,  the  computed  wall  tempera- 
ture becomes  the  temperature  of  the  fluid  film  adjacent  to  the 
wall  and  the  coolant  temperature  becomes  the  coolant  discharge 
temperature.  Therefore,  from  equation  67,  the  fluid  film  tem- 
perature is 


where  f(v)  is  defined  by  equation  62.  The  variable  T is  the 
fluid  temperature  adjacent  to  the  cooling  film  as  dem &£- 
strated  by  the  temperature  profiles  shown  on  figure  2 (b). 


In  evaluating  the  variables  of  equation  75,  the  characteristics 
of  the  plate-fin  material  configuration  selected  for  the  heat  ex- 
changer must  bo  known.  The  characteristics  of  several  types  of  con- 
figurations are  given  in  references  5 and  6.  Heat  transfer  characteris- 
tics are  presented  in  the  form  of  (h/GC  )(Pr)2/3  versus. R . The  nhy- 
sical  characteristics,  such  as  heat  transfer  and  flow  area,  are  also 
provided.  The  heat  transfer  data  presented  is  for  the  same  rate  of 
heat  transfer  occurring  on  both  walls  of  the  heat  exchanger.  For 
the  convection-film  cooling  method,  the  bulk  of  the  heat  flow  occurs 
across  only  one  wall.  Analytically  predicted,  the  effect  of  heat 
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flow  from  only  one  wall  is  a plate- fin  material  heat  transfer  area 
reduction  of  ‘jO  percent.  Predicted  wall  temperatures  compared  with 
empirical  data  shows  this  reduction  to  be  more  nearly  60  percent. 
Therefore, 


^Ax  * ^ ‘ ( P/aAz  -Ah  snafer/a/  hast  -tnsnr/er-  tsnra)  (8l) 


Reference  7 shows  that,  under  fixed  mainstream  conditions,  the  co- 
efficient of  convection  along  a surface  can  be  assumed  to  remain 
constant  with  or  without  a coolant  film,  i.e. 

hf  * h3  . (82) 

A coolant  flow  balance,  pressure-loss  analysis  is  conducted 
by  subprogram  CONFIM  similar  to  that  for  film  cooling  in  section 
3. 3*2.1.  The  one  difference  between  these  analyses  ia  the  pres- 
sure losses  which  occur  between  the  coolant  plenum  and  the  slot 
exit.  For  any  reasonable  design,  entrance  losses  should  be 
in  comparison  to  the  losses  which  occur  within  the  coolant  flow 
passage.  Plate-fin  heat  exchanger  characteristic  data  include 
the  Fanning  frictional  factor  (f_)  versus  R . The  flow  through 
this  system,  related  to  a series'pressure  lBss,  is 


rn&  • {(p* -*>)/[ *^ Xu  3 

+ 1 (£  )AZ  3 ft,  )]/  . (83) 


The  coefficient  (Kq(j)  Includes  all  the  pressure  losses  which  occur 
. from  the  downstream  end  of  the  heat  exchanger  passage  to  the  coolant 
discharge  slot. 
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The  subprogram  C0NFI2-1  can  either  compute  the  total  coolant 
flow  rate  or  it  can  be  a fixed  input  to  the  program.  If  the  total 
flow  rate  is  fixed,  equation  69  and  70  are  solved  to  obtain  the 
individual  slot  flow  rates.  Equations  69  through  72  are  used  to 
compute  the  overall  coolant  flow  balance  in  a manner  similar  to 
that  described  in  section  3. 3. 2.1.  To  initiate  the  computational 
procedure,  initital  estimates  are  made  of  the  heat  exchanger  fluid 
density  and  the  slot  discharge  temperature.  Initial  values  of  flow 
rate  are  computed.  The  gas  film  temperature  distribution  is  com- 
puted utilizing  equation  80.  In  turn,  using  equations  78  and  79> 
the  temperature  distribution  of  the  heat  exchanger  coolant  is  com- 
puted based  on  the  initial  estimate  of  the  coolant  discharge  tem- 
perature. The  computed  heat  exchanger  entrance  coolant  temperature 
is  then  compared  with  the  coolant  plenum  temperature  computed  from 
equation  72.  With  a mismatch  in  temperature,  an  iteration  is  car- 
ried out,  utilizing  equations  78  through  30  for  each  surface  by 
revising  the  estimate  of  coolant  discharge  temperature  until  the 
coolant  entrance  temperature  converges.  Once  completed,  the 
initital  estimates  of  heat  exchanger  fluid  density  and  slot  dis- 
charge temperature  are  revised  by  the  computed  values  and  the  com- 
plete computation  procedure  repeated,  iterating  until  temperature 
values  converge.  The  wall  temperature  distribution  for  each  sur- 
face is  then  computed  from  equation  7^  and  an  area-weighted-average 
surface  temperature  is  computed  for  each  surface  node  involved. 


The  subprogram  CONFIM  can  also  handle  the  parallel  configura- 


fv  \ 


the  exception  that  the  sign  of  the  exponential  term  of  equation  78 
is  changed.  The  computational  procedure  is  carried  out  in  a manner 
similar  to  that  for  the  counterflow  configuration  but  the  iteration 
on  coolant  slot  discharge  temperature  is  not  required. 


3. 3.2. 3 Transpiration  Cooling 

The  cooling  of  a surface  utilizing  a transpiring  gas  is  com- 
puted by  the  subprogram  TRANCL.  The  empirical  correlation  for  air 
injected  through  a porous  wall  which  was  used  in  this  program  is 
presented  in  reference  8. 

The  nomenclature  used  within  this  section  is  as  follows: 

2 

A - area,  ft. 

C - specific  heat  at  constant  pressure,  Btu/lb  . °R 

p Dl 

F - injection  mass  flow  ratio,  ( P V)/(  p V)  , dimensionless 

h - convective  heat  transfer  coefficient,  Btu/hr.ft.  R 
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- number  of  surface  increments 

- porous  wall  loss  factor,  lbf /ft.2(hr.ft.2/li  )X/n 

- mass  flow  rate,  lb  /hr. 

m 

“ porosity,  ratio  of  flow  area  to  surface  area 

- heat  flow  rate,  Btu/hr. 

- gas  constant,  lb  .ft. /lb  °R 

X LI 

- Stanton  number  (h^/ p VC^) , nondimensional. 

- temperature  or  static  temperature,  CR 

- recovery  temperature,  °R 


V 
z 

SUBSCRIPTS 

c 


o 

w 


total  temperature,  °R 

velocity,  ft. /hr. 

dimensionless  parameter 

density,  lb  /ft. 3 
m 

density  at  standard  temperature  and  pressure,  lb  /ft.^ 

m 


- coolant 

• coolant  discharge 

- malnnfi’aam  rro  c 

_ ’•  * ' 1 D""' 

“ without  transpiration  cooling 

- wall 


A sectioned  transpiring  porous  wall  model  is  shown  in  figure  4 
and  shall  be  used  for  discussion  purposes. 

The  heat  flux  by  convection  from  a gas  to  an  uncooled  wall  is 
determined  by  the  equation 


ho  (Tr  - 


(84) 


At  steady  state  conditions,  heat  added  to  the  porous  wall  from  the 
mainstream  is  removed  by  the  coolant  flow.  The  lateral  heat  trans- 
fer  by  conduction  is  assumed  to  be  zero,  i.e.  no  temperature  gradient 
ess  along  the  wall.  The  heat  addition  from  the  mainstream  is  trans- 
ferred to  and  carried  off  by  the  coolant  gas.  This  is  presented  by  the 
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heat  balance 

9 = ^AusfTr-Tu,)  = ric  Cp^  . (85) 

For  a good  design,  the  heat  transfer  efficiency  between  the  porous 
wall  and  coolant  parsing  through  the  wall  will  permit  the  coolant 
discharge  temperature  (T ) to  approach  the  wall’s  outer  surface 
temperature  (Tw) . 


Mainstream 
Flow  at 


Temperature  T 


Heat  Flux  (q/A) 

Surface  Temperature  at  Tw 

Coolant  Discharge  at  T.  / 

tc 


■Porous  Wall  Section 


Control  Boundary 

Coolant  Supply  (m  ) at  Temperature  T 
c tc 


FIGURE  4:  Transpiration  Cooling,  Sectioned  Porous  Wall  Model 


Assuming  this  to  be  the  case,  then  equation  85  can  be  written  as 

hAu(Tr-T»)  - ^Cp^(T^-To)  . (86) 


The  coefficient  of  convection  (h)  on  a transpiration  surface  is 
difficult  to  obtain.  An  equation  derived  by  reference  8,  compared 
with  the  empirical  data  presented  in  reference  9,  relates  the  heat 
flow  rate  to  a surface  with  and  without  transpiration  cooling  under 
the  same  imposed  mainstream  conditions.  Correlation  is  provided  bv 
the  relationship 


j.£_ 

z st0 


(87) 


where 

F = (a  & )/( f,  ) . 

With  the  aid  of  the  continuity  equation 
number, 


(88) 

and  the  definition  of  Stanton 


St°  = k /( Cn  f5  V.  J 

the  ratio 


(89) 


r/st.  = 


(90) 
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From  equations  84,  86,  87  and  90,  an  equation  defining  wall  ten- 
peraturc  in  terms  of  known  variables  can  be  obtained,  i.e. 

T*o  • (Tic.  •+  * 7r  ) / fz+l ) 


(91) 


where 


[($£)*(&:)']*-  Mr  . <*> 


The  subprogram  TRANCL  solves  equation  91,  with  equation  92,  for 
the  wall  temperature  distribution  based  on  the  coolant  flow  dis- 
tribution along  the  transpiring  surface. 

The  coolant  flow  conditions  are  established  by  this  program 
in  a manner  similar  to  that  discussed  in  section  3. 3.2.1.  The  same 
options  are  available,  i.e.  the  coolant  flow  rate  fixed  or  computed 
based  on  its  source  conditions.  The  pressure  loss  across  the  wall 
is  computed  from  equation  70.  Since  the  mainstream  static  pres- 
sure  may  vary  along  the  surface  being  cooled,  surface  increments 
are  taken  by  the  program  which  permit  the  coolant  flow  rate  to 
vmry  *inn2  +**  The  fW  relationship  for  the  surface 

increments,  assuming  that  the  same  plenum  feeds  the  complete  sur- 
face, is  found  to  be  * 


rf>Ci  * Aci  [&s  (Ptc-2,)  /( Be  TicJ 


Zn 


(93) 


where  the  total  coolant  flow  rate  is 


• 1 

u 


(94) 


If  the  coolant  flow  rate  is  fixed  by  the  user,  equation  93  and  94 
are  solved  to  provide  the  flow  for  each  surface  increment.  When 
the  program  is  required  to  compute  the  coolant  flow  rate,  these 
equations  in  conjunction  with  equation  71  and  72  are  solved  in  a 
manner  similar  to  that  presented  in  section  3. 3,2.1. 
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3.3*3  Exhaust  System  Radiation 

Radiant  energy  is  interchanged  between  the  internal  surfaces 
of  the  exhaust  system  and  emitted  from  the  exhaust  system  sur- 
faces to  a remote  detector.  Internal  radiation  is  required  for 
an  internal  surface  heat  balance  which  produces  the  surface  tem- 
perature distribution  within  the  exhaust  system.  The  irradiation 
incident  upon  the  detector  provides  the  signature  of  the  particular 
exhaust  system  under  investigation.  Two  parameters  offer  difficulty 
in  the  computing  of  radiation;  geometric  view  factors  and  gray  body 
interchange  factors.  This  section  presents  the  methods  by  which 
SIGNIR  calculates  these  parameters. 

The  nomenclature  used  within  this  section  and  its  subsections 
is  as  follows: 

- surface  area,  ft.2 

- matrix  element 

2 

- radiosity,  Btu/hr.  ft. 

o 

- differential  element  of  area,  ft. 

2 

- black  body  emissive  power,  Btu/hr. ft. 

- geometric  radiation  view  factor 

wwm*  4 « + n>*«Vonrro  fo 

6*  wi;  ** - .VI 

2 

- irradiation,  Btu/hr. ft. 

- heat  rate,  Btu/hr. 

- radial  distance,  ft. 

- Kronecker  delta  (1,  when  i=j ; 0,  when  i/j) 

- emissivity 

- An  angle  measured  from  the  normal  from  a surface 

- reflectivity 

SUBSCRIPTS: 

i - surface  designator  (or  matrix  row  index) 

j - surface  designator  (or  matrix  column  index) 

1 - surface  designator 

- surface  designator 


o 
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The  computing  of  net  radiation  exchange  occurring  between 
two  isothermal  surfaces,  1 and  2,  is  based  on  the  Stefan-Boltmann 
lav  sad  defined  by  the  equation 


<?;.»  = Sh  (tS-tJ)  . 


The  net  amount  of  radiaat  energy  entering  a surface  (l)  is 
the  simulation  of  the  contribution  from  all  other  surfaces,  i.e. 


<?/  * Z Q'-j  * £ s u(T>4-n*)  ■ 

J j 

Calculation  of  the  internal  geometric  view  factors  is  accom- 
plished by  subprogram  VIEW  and  the  external  view  factors  by  sub- 
program REDII.  The  systems  gray  body  interchange  factors  are  com- 
puted by  subprograms  SCRPTF  for  internal  radiation  interchange  and 
HEDI  for  the  external  case. 

3*3«3*1  Geometric  Radiation  View  Factor 


The  view  factor  between  a black  radiating  surface  (l)  and  a 
black  receiving  surface  (2)  is  that  fraction  of  the  radiation 
leaving  surface  1 that  is  Intercepted  by  surface  2.  The  defining 
equation  for  geometric  view  factors,  purely  a geometric  relation- 
ship, is  given  by 

ir  _ CoS#  COSGx  ^ (97) 

A e/A-oH*  " Fr*  ’ 


The  relationship  between  the  variables  used  and  the  surfaces  are 
shewn  on  figure  5.  The  view  factor  between  surfaces  is  obtained 
by  integrating  erver  the  surfaces,  i.e. 


rhl  „ J-  ff  C0S&1  dA, 

Aa. 


r*nst ga.  ..-  i v , vm 
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The  black  body  view  factors  are  utilized  for  gray  surfaces  when 
they  are  considered  to  radiate  in  a diffused  manner. 

The  subprogram  VIEW  computes  the  internal  exhaust  system 
geometric  view  factors  required  for  the  radiation  calculation. 

This  program  was  developed  for  use  under  NASC  contract  number 
N00019-68-C-0571  and  is  completely  documented  by  reference  10. 

The  exhaust  system  is  sub-divided  into  isothermal  segments  (nodes) 
for  the  mathematical  model.  Each  segment  is  represented  by  a 
frustum  of  a cone,  a cylinder  or  a disk  of  common  axis.  The  geo- 
metric considerations  for  these  axisymmetric  segments  substantially 
reduce  the  integration  problems  confronted  in  the  solution  of  equa- 
tion 98*  The  subprogram  VIEW  takes  advantage  of  the  system  axi- 
syimnetric  characteristic  and  can  determine  the  view  factor  between 
any  of  the  models  selected  segments.  It  also  has  the  ability  to 
determine  view  factors  where  complete  or  partial  shadowing  by  an 
intervening  surface  exists. 

The  symmetry  of  the  exhaust  system  permits  equation  97  to  be 
placed  in  cylindrical  coordinates  which  requires  integration  in 
four  variables.  In  subprogram  VIEW,  one  of  the  integrations  is 
accomplished  leaving  only  three  variables  requiring  numerical 
integration  methods.  All  geometric  considerations  and  numerical 
integration  methods  utilized  have  been  presented  in  reference  10. 

A modification  was  incorporated  in  the  interest  of  reducing 
program  compatational  time  requirements.  The  increment  size 
selection  for  numerical  integration  was  permitted  to  vary  in  a 
systematic  manner.  The  program  first  defines  the  increment  size 
as  the  node  itself  and  a value  for  the  view  factor  is  computed.  The 
node  increment  size  is  then  reduced  by  one  half  by  the  program  and 
a second  view  factor  computed.  Its  magnitude  is  then  compared  with 
the  first  value  obtained  and  if  the  percent  change  in  magnitude  is 
found  greater  than  0.1%,  the  program  again  reduces  the  Increment 
size  by  one  half.  A value  for  the  view  factor  is  computed  and  a 
comparison  made  with  the  last  value  computed.  This  process  is 
continued  until  the  tolerance  requirement  is  met  or  until  the 
view  factor  is  computed  based  on  a maximum  of  16  node  increments. 
This  modification  proved  to  substantially  reduce  the  required  com- 
puter time  with  effectively  no  reduction  in  accuracy. 

3. 3*3.2  External  Geometric  Radiation  View  Factors 

The  view  factors  between  .the  detector  in  space  and  the  exposed 
exhaust  system  surfaces  are  computed  by  the  subprogram  REDII;  a pro- 
gram developed  to  support  the  analytical  computations  accomplished 
under  NASC  Contract  N00019-68-C-0571.  The  energy  collection  axea  of 
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the  detector  is  considered  to  be  am  element  of  one  square  unit 
for  calculation  purposes.  With  this  single  consideration,  equa- 
tion 97  can  be  stated  as 


-At 


<3as&,  <2ois&e. 

nr* 


JAl 


(99) 


The  normal  to  the  detector  area  is  assumed  to  intersect  the  exhaust 
system  at  the  point  where  the  exit  plarne  and  system  axis  intersect. 
Each  node  within  the  exhaust  system  is  subdivided  into  elements  such 
that  the  summation  of  the  detector  to  element  view  factors  provide 
the  numerical  integration  for  the  view  factor  from  detector  to  node. 

Subprogram  REDII  checks  positioning  of  each  element  to  determine 
if  a view  is  possible  ( 0 < 90°)  and  checks  possible  shadowing 

by  any  Intervening  surfaces.  A discussion  of  the  geometric  considera- 
tion undertaken  by  the  subprogram  REDII  is  covered  in  references  11 
and  12. 


3*3*3. 3 Grey  Body  Interchange  Factors 

Radiant  energy  is  reflected  from  non-black  surfaces.  In  an 
enclosure,  the  reflection  of  energy  can  occur  numerous  times.  These 
multiple  reflections  of  radiation  are  accounted  for  through  gray 
body  interchange  factor  for  surfaces  which  reflect  radiation  in  a 
diffused  manner.  These  factors  are  computed  within  subprogram  SCRPTF 
for  the  internal  exhaust  system  and  by  subprogram  REDII  for  system  to 
detector  irradiation. 

The  radiosity  of  an  Isothermal  surface  i is  equal  to  the  sum  of 
the  radiation  emitted,  reflected  and  transmitted.  Opaque  surfaces 
such  as  those  used  within  exhaust  systems  transmit  no  radiation, 
therefore 

‘ Ai  y £bi  (100) 


The  total  irradiation  on  surface  i is  the  sum  of  the  contributions 
from  each  of  the  other  surfaces,  i.e. 
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which  involves  the  diffuse  geometric  radiation  view  factor. 
Written  in  matrix  notation  (i  = row  index,  j = column  index), 
equations  100  and  101  become 

{Blj  ' {/  '-  fa  * fat£fa 


{fa  - [Ffafafa 


Combining  equations  102  and  103  and  solving  in  terms  of  the 
radiosity,  i.e. 

{fa  ' fafalfa  * fa£fa  , 
fa  ■ & fa  { fa  * fa  Ffa 

and 

fa}  ■ [Sij  -i°i  fa'fai  fa  % 

the  radiosity  can  be  expressed  as 

{fa  - faj  fa  fafa 
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where 


[Zij]  ■ [g{j  -ftF-.;]"  ' 


(108) 


The  net  rate  at  which  radiation  enter*  - 
area  and  time  is  equal  to  the  dif?erJn.  SY  SUrface  Per  unit 
and  the  emitted  radiation.  Assuming  thf  b®tween  the  absorbed 
surface  equals  its  emissivity  thi ^€  ^e  absorptlvlty  of  the 

viTjy,  tftls  car*  be  expressed  as 


( Ac  } ' H } . 


(109) 


Combining  equations  103,  107  with  109  gives 

<?c 


f^ij}  e fcj][Acj  *;][£*}  r {fuJ  m 


(110) 


This  ran  he  written 


f<?c}  ■ fa  A{  (v,j  - £ij )]  fa.} 


(in) 


where 
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The  gray  body  interchange  factor  from  surface  j to  surface  i is 
defined  such  that  the  net  heat  exchanged  between  the  two  surfaces 
is  given  by 


Gj-i  = AJ  ' £*l  ) 


(113) 


Qj-l  ' A (£bj  ' £bL  ) 


(11U) 


The  net  amount  of  heat  entering  surface  i is  obtained  by  the 
summation  of  the  contributions  of  each  of  the  surfaces,  i.e. 


Ql  * Y,  Q/i  ' X Ai  tfj  ~ £bz ) . 


(H5) 


Transforming  to  matrix  notation,  equation  115  becomes 


f Q } * [ A ] f^bij  - Ai  JFt>£  J (116) 


or, ' can  be  written  as 


{<?i}  ■ [a  4 - %j  £ a t*it!  ] {b„.  } . (W) 
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Comparing  the  elements  of  equations  111  and  117,  it  may  be 
•seen  that  for  i { j, 


A- 


= A A 


(no) 


and  for  i * j, 

A r £ +6;Ac(T>U~l)  , (119) 


Since  each  of  the  surfaces  considered  is  isothermal,  the  inter- 
change factor  of  a surface  from  itself  is  not  required  and  the 
diagonal,  equation  119  need  not  be  solved.  Equations  118,  along 
with  108  and  112,  are  solved  to  obtain  the  gray  body  Interchange 
factors  required  for  internal  radiation  calculations. 

The  interchange  factors  for  the  detector  involve  only  the 
irradiation  from  the  exhaust  system,  consequently  the  Kronecker 
delta  terms  of  equations  111  and  117  are  zero  and  again,  only 
equation  118  is  reauired. 

3.3.^  Exhaust  System  Heat  Balance 

A thermal  heat  balance  is  conducted  within  the  exhaust  sys- 
tem to  provide  the  temperature  distribution  along  all  exhaust 
system  surfaces.  Each  surface  has  been  subdivided  into  isothermal 
segments  (nodes)  and  it  is  the  temperature  of  these  segments  that 
•;ia  desired.  The  subprogram  TAS1B  is  employed  to  compute  these  tem- 
peratures. This  subprogram  is  a modified  version  of  the  existing 
program  Thermal  Analysis  System  1 written  by  the  Jet  Propulsion 
Laboratory,  California  Institute  of  Technology  for  NASA  under  con- 
tract no.  NAS  7-100. 

The  nomenclature  used  within  this  section  is  as  follows: 

2 

A - surface  area,  ft. 

C - overall  conductance  coefficient,  Btu/hr.  ft.  °R 

& - gray  body  radiation  interchange  factor,  nondimension al 


'm*m&m**&mm*mmrn-i  nim*mmmi  , , «^*wr -i<i«rrii,_ 
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R - vector 

S - square  matrix 

T - temperature,  °R 

q-  - Stefan -Boltzmann  constant. 


0.1713  x 10"8  Btu/ft.2 


hr.  ( r) 


For  any  given  internal  exhaust  system  node,  a steady  state  energy 
balance  can  be  stated  as  the  sum  of  the  energy  addition  is  equal 
to  the  sum  of  the  departing  energy.  The  balance  for  a node  i can 
be  written  as 


£ C£J  A { (Ti-Tj  ) 4 £ <5Ai%.  (Ti4-  Tj4)  - O 

J J 


(120) 


/ 

I 

P« 


and  stated  as  the  net  energy  exchange  from  node  i is  zero. 

For  a system  of  nodes,  the  energy  balance  can  be  written 
in  matrix  form  as 


[cisAi]{T£}  4 [°A%j]fc*}  -o  . 


The  solution  used  follows  the  Newton-Raphson  iteration  method 
for  nonlinear  algebraic  equations  utilizing  the  matrix  form 


fadfai  'fe‘  J 


(121) 


(122) 


The  off-diagonal  elements  of  the  square  matrix  are 


*</  “ 


CijAi  -<■ 


(123) 
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while  diagonal  elements  are 


S{{  * Sii  - CijAi  ~ 4c  AL  £tCJ  Tj  , <«*>> 


fhA  vector  R has  the  form 


- 3(5  Ac  (tc4 


(125) 


Subprogram  TAS1B  solves  this  matrix,  equation  122,  for  the  tempera- 
ture (T)  vector  in  an  iterative  manner  revising  the  values  of 
temperature  in  matrix  S and  R with  each  iterative  step  until  the 
process  converges  on  a final  solution. 

3.3.5  Exhaust  System  Radiation  Emission 

The  spectral  radiation  passed  from  an  engine  exhaust  system 
to  a detector  is  dependent  on  the  conditions  of  the  systems  exposed 
hot  metal,  the  separation  between  the  two  and  the  media  through  which 
the  energy  must  pass.  The  following  section  presents  the  theory 
utilised  in  computing  the  spectral  emission  from  an  exhaust  system 
to  a detector  and  the  attenuation  of  infrared  energy  as  it  passes 
through  the  atmosphere. 

3. 3. 5,1  Spectral  Emission 

- Th*  spectral  radiation  passing  from  the  exhaust  system  surfaces 
to  a remotely  located  infrared  detector  is  computed  by  the  subpro- 
gram REDI.  This  program  was  developed  by  VAC  under  NASC  Contract  No. 
S00019-68-C-0571  and  reported  in  reference  11.  The  only  modification 
incorporated  by  this  program  is  a reduction  in  the  size  of  the  spectral 
wavelength  band  interval  resolved. 
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,.  . The  following  is  a list  of  nomenclature  used  throughout 
this  section. 


A 

B 

C 

exp 

■b 

F 

& 

Q 

r 

R 

T 

U 

Yfr) 

Sij 

A 

€■ 

A 

P 

fl 


- surface  area,  ft.2 

- radiosity,  Btu/hr.  ft.2 

- a constant 

- exponential  (e) 

- black  body  emissive  power,  Btu/hr.  ft.2 

- geometric  radiation  view  factor 

- gray  body  interchange  factor 

- heat  rate,  Btu/hr. 

- range,  ft. 

- a constant 

- temperature,  °R 

- a constant 

- Gauss  quadrature  formula  coordinant 

- Kronecker  delta  ( * 1,  when  i = j;  =0,  when  i / j) 

- incremental  change  in  a variable 

- emissivity 

- wavelength,  microns 

- reflectivity 

- solid  angle,  steradians 


SUBSCRIPTS: 

i - matrix  row  index  (or  surface  indicator) 
j - matrix  column  index 


The  amount  of  radiant  energy  in  the  wavelength  band  fran  ^ ttf 
** **  which  is  incident  on  an  IR  detector  may  be  written  as  the  sum 
of  the  energy  contributions  from  each  of  the  individual  exhaust  sys- 
tem surface  nodes.  That  is, 


(-) 

1 A '?  avAa 


* I n 


A -6a  A-/ A A 


(126) 
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The  radiosity  may  be  expressed  in  terms  of  the  black  body  emissive 
power  for  each  surface.  The  expression  is  obtained  in  a manner 
similar  to  that  of  equation  106.  Written  in  terms  of  matrix  no- 
tation, 

{&i4  a xja*}  m [&j  ftjj  ^ 'bij  a *•  a+A>  } , 


where  the  radiosities  are  shown  to  be  independent  of  the  properties 
of  the  detector.  Equation  126  and  127  are  combined  to  form  the 
matrix  equation 


/(f>L »**>}  * 


From  a development  similar  to  that  used  in  section  3«3«3«3»  equa- 
tion 128  can  be  written  as 

{(fi {*}}{*«.»  + *“*} 


where  is  the  gray  body  interchange  factor  from  the  receiver  to 
an  exhaust  system  surface  node. 

The  black  body  emissive  power  of  surface  node  i is  obtained  by 
integrating  Planck's  equation  for  radiation  between  the  limits  of 
a to  phap  where  the  value  of  AP  * 0.05  microns.  That  is, 


Pi  -Hf  Pi+APi 


C,  d* 

?s  [*xp(Cz  /pTc  ) " tj 


(130) 
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The  integral  is  evaluated  by  means  of  numerical  integration  with 
the  Gauss  quadrature  formula  for  four  selected  abscissas.  The 
following  equations  are  employed. 


/ e,  y, Y?, * £*  Yt'M)]  (131) 


where 


a, 


and  = 

U2  = 

U3  = 
U4  “ 


u [ 

<zxp  (e,  ) - / J , 

(132) 

* 71  V (Vn 

■v  •/*)(&*) 

(133) 

-0.4305682 

= 0.1739274 

(134) 

-0.1699905 

Rg  = 0.3260726 

0.1699905 

r3  = 0.3260726 

0.4305682 

r4  = 0.1739274 

The  subprogram  REDI  computes  the  spectral  radiation  incident 
on  the  detector  by  solving  the  matrix  equation  129  utilizes 
equation  131  to  compute  the  black  body  spectral  emissive  power. 
This  radiant  intensity,  incident  radiation  per  unit  solid  angle, 
is  computed  by  multiplying  the  spectral  radiation  flux  of  equa- 
tion 129  by  the  square  of  the  range.  The  solid  angle  itself 
would  be  calculated  as 


1 


fl  r A / rl 


(135) 


The  computed  radiant  intensity  is  therefore  independent  of  the 
energy  collection  area  of  any  specific  detector. 
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3A  SPECIAL  COMPUTATION  CONSIDERATIONS 

. ..  ft'ogram  SIGNIR  is  basically  self  contained,  dependent  only 
on  the  system  model  produced  from  the  physical  exhaust  system  con- 
figu,  ation.  With  the  latitude  provided  in  exhaust  system  configura- 
tions, certain  specific  items  concerning  a system  may  be  required  to 
be  provided  by  the  user  to  complete  the  model  and  derive  the  desired 
information.  The  following  sections  are  included  to  assist  the  user. 
They  cover  a method  to  approximate • system  thrust  change  and  establish 
the  required  form  for  pressure  loss  and  heat  transfer  parameters 
necessary  for  program  input. 
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3.4.1  Thrust  Change  Resulting  From  Suppressing  Exhau3t  System 

A method  is  presented  by  which  the  exhaust  system  thrust 
change  associated  with  the  infrared  suppression  of  a system  can 
be  approximated  by  the  user.  To  be  capable  of  obtaining  this  Infor- 
mation, both  suppressed  and  unsuppressed  configurations  must  be  run 
by  program  SIGNXR.  Available  as  output  data  from  each  of  these  con- 
figurations is  a computed  ’’thrust  loss  factor."  This  factor  is 
the  summation  of  the  internal  exhaust  system  axial  pressure-area 
forces,  surface  friction  forces  and  the  coolant  momentum  flux 
gain.  The  method  assumes  that  the  suppression  of  the  system  does 
not  effect  engine  operating  conditions  and  the  flow  condition 
entering  the  exhaust  system  are  unchanged  by  the  configuration 
change. 


The  nomenclature  used  within  this  section  is  as  follows: 


A 


M 

e 

m 

c 


n 

n 

P 

q 

r 

S 

(3 

TT 

u 

X 


cross  sectional  area,  ft. 
coefficient  of  discharge 
differential  of  a quantity 
thrust  force  factor,  lb. 
force  vector,  lb. 
coolant  duct  drag  force,  lb. 
coolant  duct  reaction  force,  lb. 
frictional  force,  lb. 
pressure-area  force,  lb. 

Mach  number 

coolant  mass  flow,  lb/sec. 
unit  normal  vector 
number  of  coolant  ducts 
pressure,  paia 
dynamic  pressure,  psf 
radius,  ft. 

2 

projected  coolant  duct  area  normal  to  flow,  ft. 

reaction  force,  lb. 

velocity  vector,  ft. /sec. 

velocity,  ft. /sec. 

system  coordinate 


Cr 


51 





ww* 


p 


V). 


VOUGHT  AERONAUTICS  COMPANY 


- ratio  of  specific  heats 

- change  in  a parameter 

3 


r 

A 

/o  - density,  lb/ft. 
SUBSCRIPTS 


amb 

cd 

ce 

cs 

ent 

ex 

sup 

tc 

unsup 

Special  notations: 


ambient 

coolant  system  discharge 
coolant  system  entrance 
control  surface 
entrance 
exit 

suppressed 

tailcone 

unsuppressed 


applies  to  tailpipe 
applies  to  centerbody 
applies  to  splitter 


The  thrust  acting  on  a system  is  determined  from  the 
momentum  equation  for  steady  state  flow  conditions,  i.e. 


u (/ou-ri)  c/A  - 


(159) 


Figure  6 shows  components  of  an  exhaust  system  similar  to 
the  demonstration  problem  presented  in  Appendix  III.  The  components 
have  been  isolated  such  that  individual  force  balances  can  be 
performed  on  each  of  the  surfaces.  A force  balance  conducted 
in  the  axial  direction  on  the  exhaust  system  tailpipe  (Figure 
6a)  produces 


■/  ( ■+  (fp*)*  ' (A^i  - A**) 


(166) 
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' ' ' 1 lMrw’ 


.here  the  forced'  Is  the  reaetlon  force.  A force  balance 
Conducted  on  the  exhaust  system  centerbody  (figure  6b) 
suits  In  the  equation  form 

£r  + Ov")x  *(?£),  ‘1  rifled' 

-In  addition  to  these,  a force  balance  on  the  splitter  with  the 
coolant  ducts  (figure  6c)  yields  the  equation 

+ ( FJ" J,  ■+  (*%“)*  * (Fdr^  )ed  -(p°)x. 

= - (ri±u)c<e 

The  total  thrust  acting  on  the  system  is  equal  to  the  sum 
of  these  individual  forces  (thrust  components), 

If" 

With  equations  159  through  163,  the  total  system  thrust  becomes 

' - 1 (f=/K  -Z(F»)*  + P*r*b(A*"+~A**) 

- (PA)ic  + [l(rhQU)cd  - 
- C m<>  U h ] 


* 
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If  the  engine  dimensions,  ambient  pressure  and  the 
pressure  applied  to  the  tailcone  are  assumed  to  remain 
fixed  for  configuration  comparison,  the  terms  Pftrnh  ana  V 

are  the  same  for  the  two  configurations.  The  program  computes 
a "force  factor"  which  is  equated  by  the  expression 

p =•  -z  (Ff>*  - Z (?/»>*  * £ (^u)od ' 


(165) 


The  thrust  change  between  configuration  can  be  written  as 

A £/  ~ ^unsup  ~ ^sup 


(166) 


and  expanded  to  the  form 

Atf  ’ [P  - ' [F  ' 


sup 


-b 


+ ( ft1  c U )CtT  Jsup 


This  equation,  developed  from  a rather  specific  configuration, 
is  general  in  nature  and  should  be  found  to  apply  to  almost  all 
jet6 configurations  encountered.  In  cases  where  configur  - 

tions  havethe  same  tailpipe  exit  area  and  no  cooling  system 
exists  on  the  suppressed  configuration,  equation  1 7 re  uce 
the  simple  form  of 


A - C^hnsup  ~ (P*  ^ 


'sup 


(167) 


(168) 


For  this  particular  case,  the  thrust  change  is  only  dependent  on 
the  "force  factors"  provided  by  the  program  SIGNIR. 
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In  genera,  the  user  will  have  to  solve  the  thrust  change 
equation  form  shown  in  equation  167.  With  the  program  input  and 
output  data,  the  first  two  terms  of  this  equation  can  be  evalua- 
ted. The  last  term  involves  a knowledge  of  the  particular  sup- 
pressed configuration  under  investigation.  If  coolant  ducts  are 
present  and  pass  through  the  exhaust  gas  stream,  the  drag  force 
produced  is  a function  of  the  duct  configuration  and  the  gas  stream 
properties.  This  force  must  be  computed  by  the  program  user.  To 
assist  in  this  effort,  the  parasite  drag  caused  by  ducts  can  be 
evaluated  by  the  equation. 


= C 


■p  h S Q — Cp  ft  S 


rPM z 


The  values  of  pressure,  Mach  no.  and  ratio  of  specific  heats  are 
available  from  the  program  for  the  specific  gas  stream  location 
required.  The  drag  coefficient  is  dependent  on  duct  profile  and 
flow  Reynolds  number. 

The  coolant  duct  entrance  momentum  flux,  (mu),  must  be  in- 
cluded if  coolant  air  is  extracted  from  the  system  as  shown  in  the 
sample  problem.  If  test  data  is  not  available,  this  function  can 
be  estimated  by  assuming  the  entrance  is  a fairly  well  designed 
inlet  whose  flow  contraction  ratio  is  approximately  0.7.  The 
momentum  flux  of  the  flow  stream  captured  by  the  coolant  system 
entrance,  is  approximated  by  the  equation 

(ricuL  — O.lTnPAM*'  { 


where  the  value  of  pressure,  Mach  no.  and  ratio  of  specific  heats 
are  provided  by  the  program  as  stream  conditions  at  the  coolant 
duct  entrance  location. 


The  thrust  change  associated  with  the  infrared  suppression  of 
an  exhaust  system  can  be  computed  from  equation  167.  Specific  in- 
formation required  to  accomplish  this  task  can  be  obtained  from 
equations  169  and  170  and  the  output  data  from  the  program  SIGKIR. 
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3.4.2  Pressure  Loss  Parameters 


When  surface  cooling  is  employed  by  the  exhaust  system,  the 
pressure  loss  characteristics  must  be  supplied  by  the  program  user 
to  define  the  particular  cooling  system  configuration  employed. 
Program  input  data  requires  values  of  the  loss  coefficient  (K)  and 
the  flow  exponent  (n)  for  specific  segments  of  the  cooling  system. 
If  these  values  are  not  available  from  measured  data,  they  must  be 
acquired  through  empirical  data  presented  in  available  reference 
material  such  as  reference  16. 


is 


The  nomenclature  used  within  this  section  is  as  follows: 


The 


A 

K 

K 

& 


n 

9 


u 


(S' 

form  of 


- flow  area,  ft. 

- loss  coefficient,  psf/(lb./sec. )n 

- nondimensional  loss  coefficient 

- mass  flow  rate,  lb/sec. 

- flow  exponent 

2 

- dynamic  pressure,  lb/ft. 

- flow  velocity,  ft. /sec. 

- change  in  total  pressure,  lb/ft.  . 

- fluid  density,  lb/ft.^ 

- density  at  standard  atmospheric  temperature 
and  pressure,  lb/ft.3 

- fluid  density  ratio  ( ) 

the  pressure  loss  equation  used  by  program  SIGNIR 


oAPt  * x (m)n  (171) 

in  which  K is  a dimensional  loss  parameter  and  n af  2.  These 
values  can  be  obtained  directly  from  system  test  data.  When 
flow  loss  estimations  are  required,  other  empirical  loss  data 
must  be  utilized.  These  are  frequently  presented  as  a dimension- 
less loss  parameter  (K.)  versus  the  flow  Reynolds  number.  The 
pressure  loss  associated  with  this  loss  parameter  is  defined  by 
the  equation 


APt  • , 


(172) 
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where 


? * i /°uZ  • 


(173) 


Using  a value  of  two  for  exponent  n,  these  two  loss  parameters 
are  related  by 


* ' **/  ( */*&/>*) 


(174) 


where  the  area  (A)  is  that  associated  with  the  dynamic  pressure 
( £ ) used  in  equation  172.  For  pressure  losses  which  occur  in 
series,  the  loss  coefficient  is  computed  by  the  equation 


* “ * 


(175) 


When  a parallel  flow  system  is  employed,  the  required  loss  co- 
efficient. Is  nomnuted  hv  the  emiatlon 


* 7-2 


a an 


(176) 


Equations  175  and  1?6,  based  on  n = 2,  provide  the  user  a means 
of  computing  the  loss  coefficient  for  systems  in  which  estimated 
values  are  required. 

3.4.3  Heat  Transfer  Parameters 

Special  paths  of  heat  flow  for  the  exhaust  systems  mathemati- 
cal model  caui  be  established  by  the  user  to  supplement  the  batic 
program  process.  These  improve  the  programs  general  heat  balance 
by  incorporating  the  characteristics  of  a specific  system  through 
defining  additional  heat  transfer  paths  to  heat  sources  or  sinks. 
Heat  transfer  parameters  are  used  to  define  the  barrier  to  heat 
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transfer  for  any  surface  coolant  supply  system  and  to  special 
fluid  nodes  (see  section  4.2.6).  In  both  cases,  an  overall 
heat  transfer,  coefficient  is  required  by  the  program  to  describe 
this  barrier  to  heat  flow  if  it  is  to  be  included  in  the  heat 
balance. 


The  following  is  a list  of  nomenclature  used  within  this 
section. 


A 

a 

b 

c 

h 

k 

L 

Pr 

Re 

t 

U 


area  normal  to  heat  flow  path,  ft. 
Reynolds  number  - exponent 
Prandtl  number  - exponent 
convection  coefficient  constant 


2 o 

- convection  film  coefficient,  BTU/hr.  ft.  R 


condu''  ivity,  BTU/hr.  ft.  R 


characteristic  length,  ft. 
Prandtl  number 
Reynolds  number 
material  thickness,  ft. 


2 0 

- overall  heat  transfer  coefficient,  BTU/hr.  ft.  R 


SUBSCRIPTS: 

d - duct 

i - inner 

o - outer 


For  a surface  cooling  system,  an  overall  heat  transfer  co- 
efficient (UA)  can  be  employed  to  establish  the  barrier  for  a heat 
path  from  a known  heat  source  to  the  coolant  passing  through  the 
coolant  supply  ducting  system.  If  the  coolant  flow  rate  is  pro-  ^ 

vided  as  a fixed  program  input,  a coefficient  is  not  required.  If 
the  heat  source  is  such  that  no  heat  transfer  will  occur  (source 
and  coolant  temperatures  identical),  a UA  value  of  one  will  permit 
proper  program  operation  without  requiring  the  user  to  compute  the 
actual  value.  When  a temperature  difference  exists  between  the 
coolant  and  heat  source,  a supply  system  overall  heat  transfer  co- 
efficient should  be  computed  by  the  user  for  program  input.  The 
path  from  the  heat  source  to  coolant  will  normally  pass  through  a 
series  of  barriers  which  are  both  the  inside  and  outside  gas  films 
and  duct  conduction.  For  this  series  configuration,  the  overall 
heat  transfer  coefficient  can  be  computed  by  the  equation 


-u.  W i-i,  ' 
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The  value  of  conduction  (k)  for  the  duct  is  a material  property 
and  can  be  found  in  reference  material  such  as  reference  16.  Con- 
vection film  coefficients  are  functions  of  gas  flow  properties  and 
the  duct  configuration.  These  are  normally  evaluated  through  an 
equation  of  the  form 


h ■ c f(Rtf  C8-)‘‘  _ (178) 


Values  for  the  constants  of  a,  b,  and  c axe  usually  presented  as 
a function  of  Reynolds  number  for  a specif i’  configuration  and 
can  be  found  in  numerous  references;  two  common  references  are 
17  and  18. 

When  special  fluid  nodes  have  been  incorporated  in  the  exhaust 
system  model,  an  overall  heat  transfer  coefficient  is  required  to 
define  the  barrier  to  heat  transfer  between  the  fluid  node  and  its 
adjacent  node  in  the  system.  These  normally  will  follow  a path 
from  a surface  to  fluid  node  passing  through  a conducting  wall  and 
a convection  film.  The  overall  heat  transfer  coefficient  is  com- 
puted in  manner  similar  to  that  presented,  in  this  particular  case 
the  equation  is  found  to  be 

U/Q  X ^ h~ ) . (179) 

Again,  the  values  of  k and  h must  be  obtained  based  on  empirical 
data. 

With  these  equations,  the  user  can  provide  the  overall  heat 
transfer  required  as  program  input  in  most  of  the  cases  that  will 
confront  him. 
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5.0  PLUME  SIGNATURE 

The  following  pages  are  reproductions  of  documents  which  are  difficult 
to  obtain  or  considered  essential  to  the  definition  of  the  plume  emission/ 
transmission  calculations . 
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THE  UNIVERSITY  OF  MICHIGAN 
THE  INFRARED  AND  OPTICS  LABORATORY 

MEMO  TO:  File  28  January  1970 

FROM:  G.  Lindquist 

SUBJECT:  How  to  use  FS+TAU5. 

FS  is  a general  program  to  computer  radiances,  absorptances,  and  apparent 
temperatures  through  a path  having  any  desired  temperature  and  pressure  variations 
with  distance  and  can  be  used  for  isothermal  and  nonisothermal  single  gases  or  three - 
gas  mixtures,  two  of  which  can  be  absorbing.  The  latest  band  model  techniques  are 
used  and  the  lines  are  allowed  to  be  overlapping  (although  an  option  to  use  the  non- 
overlapping  line  case  is  available).  An  option  is  also  available  which  allows  the 
radiances  from  a hot  gas  source  to  be  first  computed  and  then  adjusted  by  the  trans- 
mittance of  an  atmospheric  path  inserted  in  front  of  the  hot  gas  source.  This  method 
is,  of  course,  incorrect  and  is  used  only  for  comparison  purposes. 

The  path  through  which  the  calculation  is  to  proceed  is  described  by  5 linear 
arrays  XTEMP,  TEMP,  PRES1,  PRES2,  andPRES3.  XTEMP  describes  distance  from 
the  observer.  The  other  four  arrays  describe  the  temperature  and  the  pressures  of 
absorbing  gas  1 and  2 and  the  broadening  gas,  3,  at  the  distances  corresponding  to  the 
values  of  XTEMP  (see  Figure  1).  Note  that  the  observer  is  always  at  a value  of  XTEMP 
= 0.  NT  specifies  the  number  of  values  that  are  in  each  of  these  arrays. 

The  size  of  the  increments  through  which  the  integration  over  the  above  described 
path  takes  place  is  determined  by  the  four  variables  NDIV.DL,  PATHLN,  and  IPL.  Two 
options  are  available  determined  by  the  value  of  IPL.  For  IPL  = -1,  PATHLN  is  taken 
to  contain  the  total  distance  over  which  the  integration  is  to  take  place  (its  value  must 
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be  s the  last  value  in  the  array  XTEMP) . This  path  length  is  then  divided  into  NDIV 
equal  increments  of  length  DL(1)  = ' ’ • = DL(NDIV)  = . For  IPL  = +1, 

the  value  of  PATHLN  is  ignored  and  NDIV  values  of  the  incremental  length  DL,  must 
be  read  in.  The  temperatures  and  pressures  for  each  segment  of  length  DL  are  taken 
by  interpolation  in  the  functional  variations  of  temperature  and  pressure  described 


by  TEMP,  PRES1,  PRES2,  and  PRES3  versus  XTEMP.  In  this  program  all  distances 


o 

are  in  centimeters,  all  temperatures  in  K,  and  all  pressures  in  atmospheres. 

The  frequency  range  over  which  the  calculations  are  made  is  determined  by 
the  three  variables  NFREQ,  NUFRST,  and  NUINC.  NUFRST  describes  the  first 
frequency  to  be  computed,  NUINC  describes  the  amount  by  which  the  frequency  is 
to  be  incremented  between  successive  frequencies  for  which  computations  are  de- 
sired and  NFREQ  is  die  number  of  frequencies  for  which  computations  are  desired. 

If  NUFRST  is  set  to  be  less  than  0,  then  the  program  expects  to  find  the  NFREQ 
frequencies  for  which  calculations  are  desired  in  the  array  NU. 

The  band  model  parameters  for  the  absorbing  gas  or  gases  are  expected  to  be 
in  an  external  file  and  are  to  be  in  the  proper  order.  The  array  GAS  contains  the  names 
of  the  absorbing  gases,  and  if  GAS(1)  = 'H26'  (as  is  set  in  the  data  statement)  the  first 
file  must  contain  die  band  model  parameters  for  H20,  etc. 


The  band  model  parameters  used  are  die  average  line  width  to  spacing  ratio 
multiplied  by  ~ called  /?,  and  the  average  line  strength  to  spacing  ratio,  called 
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Line  # 

1 


2 


3. 


>1 

ii  i i*mtam-r^m*rx,. 
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Contents  of  Line 

Number  of  temperatures,  NTEMP,  for  which  p is  read  in, 
number  of  frequencies , NNU,  for  which  p is  read  in 
(format  = 215) 

Up  to  72  characters  of  format  to  be  used  to  read  the  data 
immediately  following  this  line  (format  = 18A4) 

Values  of  frequency  for  which  p is  given . There  should  be 
as  many  values  as  specified  in  line  i *. 

Up  to  72  characters  of  format  to  be  used  to  read  the  data 
immediately  following  this  line  (format  = 18A4) 

Values  of  temperature  for  which  p is  given.  There  should 
be  as  many  values  as  specified  in  line  1 . 

Up  to  72  characters  of  format  to  be  used  to  read  the  data 
immediately  following  this  line  (format  = 18A4) 


• Values  of  p as  a function  of  frequency  and  temperature . The 

order  of  the  values  is  as  follows:  values  of  p at  the  first  fre- 
quency for  the  NTEMP  values  of  temperature,  the  values  of  p 
at  the  2nd  frequency  for  the  NTEMP  values  of  temperature, 
up  through  die  NNU  frequencies . 

The  same  series  of  lines  is  used  to  read  in  the  k’s  as  a function  of  temperature  and 
frequency  (starting  with  a line  similar  to  line  1). 
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Currently  the  band  model  parameters  for  water  vapor  are  contained  in  a file 
named  SX90:BKH20  and  the  corresponding  parameters  for  C02  are  contained  in  SX9Q: 
BKC02.  The  program  values  have  been  preset  so  that,  unless  the  user  reads  in  other 
values,  GAS(l)  is  ' H20’  , GAS<2)  is  ’ C02’  and  GAS(3)  is  ' N2’  . Appropriate  values 
for  the  broadening  coefficients  are  also  preset.  The  upper  and  lower  frequency  limits 

for  these  values  are  1750  and  4300  cm’1.  They  are  valid  for  temperatures  between 
300°K  and  3000°K . 


The  equation  used  for  the  effective  half -width  is: 
n 


^e  = 


ways: 


y k(^,T.)  x (273.15/T)  X /?(y,T)  X (P.  + a x P ) X (DL) 
1=1  1 1 i B,i  ' 7i 

~n~  — — _ 

£ k(i;-Ti)  X (273.15/Tj)  XP.X  (DL). 

The  exponent  (rj)  of  the  overlap  factor  {fi/fif  , can  be  treated  in  three 


( 1)  set  equal  to  namelist  input  EXPO 

(2)  calculated  as  f(xe)/yxe  X .6366 

(3)  calculated  so  that: 

rj  = 0.  for  x s 1. 

JL 

V = (XL  ' l-)/9.  for  1.  < x < 10. 
r?  = 1.  for  xL  a 10. 

P,.  T4  and  Pgl  are  d,e  pressure  of  the  absorbing  gas,  the  temperature  and  the  pressure 
of  the  broadening  gas  respectively  corresponding  to  the  i-th  value  of  DL.  In  this  program 
GAS(3)  is  always  treated  as  the  broadening  gas  , t,  is  fte  broadening  coefficient. 
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The  input  parameters  are  read  from  two  logical  10  units:  the  band  model 
parameters  are  read  from  unit  3 and  all  other  input  data  from  unit  1 . The  data 
read  from  unit  1 are  read  in  via  a namelist  statement  whose  name  is  IN.  The  pro- 
gram object  module  is  currently  in  two  files  which  must  concatenated  to  run  the 
program . They  are  SX90:NFSOBJ  and  SX90:NTAu5OBJ.  Subroutines  from  the  file 
ST82.SUBPACK  are  also  required.  All  printed  output  is  written  on  logical  10  unit 
6.  Punched  cards  for  the  Benson -Lehner  plotter  (if  desired)  are  written  on  unit  2. 

The  logical  10  unit  on  which  data  for  the  CALCOMP  plotter  are  written  can  be 
specified  by  the  value  of  ICPLOT . Thus  a sample  run  statement  appears  as  follows 
(assuming  it  is  being  run  from  project  SX90).  $RUN  NFSOBJ+  NTAU50BJ+  ST8 2: SUBPACK ; 
1 = <Data>  2 = <Hot>  3 = BKH20+BKC02  6 = <Output> 

where  <Data  > is  the  name  of  the  file  on  which  the  users  data  for  the  namelist  IN 
is  stored 

<Plot>  is  the  name  of  the  file  where  plotter  card  images  for  the  Benson - 
Lehner  plotter  are  to  be  written.  The  2=  <Plot>  is  not  needed  if 
no  plotter  cards  are  required. 

<Output>  is  the  file  where  line  images  for  the  printer  are  to  be  written. 

Following  is  a glossary  of  the  parameters  that  can  be  varied  via  the  name- 

list  IN. 


GAS 

NT 


name(s)  of  gas(es),  format  A4 

number  of  input  temperatures  (integer)  and  pressures  to  be 


read  in . 
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TEMP 

XTEMP 

PATHLN 

NFREQ 

NUFRST 


NUINC 

NU 

NDIV 


NAG 

ALPHA 


\ 
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Input  temperatures  (°K)  (real  numbers) 

Path  positions  for  temperatures  and  pressures  (real  numbers)  (cm) 
Total  path  length  (cm)  (real  number)  needed  only  if  IPL  < 0 
Number  of  wavenumbers  for  which  the  calculations  are  to  be 
made  (integer) 

First  wavenumber  for  which  the  calculation  is  to  made  (real 
number).  If  NUFRST  £ 0.,  the  wavenumbers  (NrJ)  are  to  be 
read  in . 

Wavenumber  increment  (real  number).  This  is  read  in  only 
when  NUFRST  > 0.  and  the  NU's  are  to  be  calculated. 
Wavenumbers  (real  numbers);  to  be  read  in  only  if  NUFRST  £ 0 . 
Number  of  increments  over  which  the  summation  is  to  be  taken: 


(t'Ti’i) 


TAU(l/,i)) 


If  NDIV  is  negative,  NDIV  increments,  temperatures,  pressures 
of  gases  1, 2, 3 are  written.  If  NDIV  is  1,  or  -1,  the  case  is  iso- 
thermal and  the  first  values  of  TEMP  and  PRES1,  PRES2,  and  PRES3 
are  used  (integer), 
number  of  absorbing  gases  (integer) 

broadening  coefficients  (real  numbers)  e.g.  ^for^O-^; 

.769  for  CO  -N„.  They  must  correspond  in  order  to  absorbing 


gas  order. 
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PR  ESI 
PRES2 
PRES3 
IPL 

DL 

ITYPE 

EXPO 


IPRNT 

ICPLT 


Pressures  of  gas  1 (real  numbers),  (atmospheres) 
Pressures  of  gas  2 (real  numbers),  (atmospheres) 
Pressures  of  gas  3 (real  numbers),  (atmospheres) 

> 0,  DL’s  read  in  1 

' integer 

^ 0,  the  program  sets  DL(i)=PATHLN/NDIV 
Path  increments  (NDIV)  in  number)  (real  numbers).  If 
IPL  > 0,  DL’s  read  in. 
fl.  overlapping  lines 

2,  nonoverlapping  lines  (remember:  calculations  for  gas 
mixtures  cannot  be  done  with  ITYPE=2) 

t 

2 o.,  r\  = EXPO 
= -9 . , r)  = 0 . for  x s 1. 

L 

r?  = (x  -l.)/9.  for  1.  < xT  < 10. 
rf  = 1 . for  xL  s 10 . 

< 0.  and  t -9.,  tj  = f(x  )A/x  .6366 

C 0 


integer 


► real  number 


integer 


> s are  Printed  for  IPRNT  wavenumbers  and  NDIV 

increments 

s 0,  TAU’s  are  not  printed., 

b 

[<  1,  data  is  not  stored  for  CALPLOT. 

|>  1.  data  stored  for  CALPLOT  on  logical  IO  unit  number  ICPLT. 
In  this  case,  the  integer  value  of  ICPLT  must  be  set  equal  to  some 
^file  name  in  the  $ RUN  command. 
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MEMO  TO:  File 

FROM*  G.  Lindquist. 

SUBJECT:  Additions  and  Changes  to  FS 


17  April  1970 


Since  the  last  memo  was  written  several  changes  have  been  made  on 
FS+TAU5.  These  changes  make  it  more  flexible  and  update  the  procedure  so 
that  short  flames  through  long  paths  are  now  correctly  computed.  These 
changes  are  outlined  here . 

I)  A third  option  is  now  available  for  the  treatment  of  the  input  path  data . 

In  the  previous  version  (see  memo  of  28  January  1970)  the  path  element  parameters 
were  always  determined  by  interpolation  in  arrays  describing  the  temperature  and 
pressure  variation  along  the  path.  These  arrays  are  TEMP,  PRES1,  PRES2,  and 
PRES3  with  the  array  XTEMP  containing  the  path  distances  corresponding  the  values 
in  the  other  4 arrays . The  new  option,  used  by  setting  IPL  = 2,  allows  the  interpo- 
lation procedure  to  be  bypassed,  and  the  values  in  the  arrays  TEMP,  PRES1,  PRES2, 
and  PRES3  are  identified  directly  with  the  corresponding  values  in  the  path  element 
array  DL.  The  array  XTEMP  is  ignored  when  this  option  is  used.  For  example 
with  IPL  = 1,  the  values  of  TEMP  would  be  associated  with  the  values  of  XTEMP  as 
follows: 


TEMP(1)  = 300 
TEMF(2)  = 600 
TEMP(3)  = 900 


XTEMP(  1)  = 0 
XTEMF(2)  = 10 
XTEMF(3)  = 20 


which  would  describe  a path  with  temperture  that  increases  from  300°K  to  900°K  in 
the  first  20  cm.  A value  of  DL(1)  = 10.  would  result  in  linear  interpolation  in  the 
array  TEMP  to  obtain  the  value  of  temperature  for  that  DL . Thus  the  temperature 
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obtained  would  be  450°K.  For  the  same  input  data,  if  IPL  = 2,  the  value  of  TEMF(1) 
would  be  identified  directly  with  DL(1),  i.e.,  the  first  element,  10  cm  long,  would 
have  a temperature  of  300°K. 

II)  A change  is  made  in  the  way  atmospheric  paths  are  treated.  It  is  now 
possible  to  insert  an  atmospheric  path  in  front  of  the  hot  gas  and  calculate  the 
result  correctly  without  having  to  re-insert  the  hot  gas  portion  of  the  plume  along 
with  it.  To  do  this  involves  a change  in  the  use  of  the  variable  METHOD.  Previ- 
ously setting  METHOD  = ' WRONG'  set  the  program  so  that  it  would  treat 
the  input  namelists  in  pairs  and  use  the  second  namelist  to  compute  attenuation 
values  which  were  then  applied  to  the  results  from  the  first  namelist  (using  "wrong" 
method).  Setting  METHOD  = ' CORRECT'  returned  the  program  to  its  original 


mode. 


In  the  revised  program,  the  specification  of  the  variable  METHOD  either 
’ WRONG'  or  ' CORRECT'  in  a namelist  implies  that  this  namelist  is  to  be  used 
as  an  atmospheric  path  in  front  of  the  path  described  by  the  last  namelist  in  which 
METHOD  was  not  specified.  The  name  to  which  METHOD  is  set  determines  the 
type  of  calculation  to  be  done,  correct  or  wrong. 

Thus  if  the  radiance  from  a given  hot  gas  is  to  be  calculated  as  seen  through 
several  paths,  using  both  the  "CORRECT"  and  the  "WRONG"  methods,  the  name- 


lists would  appear  as  follows: 
hot  gas 


& IN  (Description  of  hot  gas)  & END 


Correct  adjustment  of  4 IN  METHOD  = ' CORRECT'  , (Description  of 

hot  gas  by  first  path:  first  path)  4 END 
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Correct  adjustment  of  hot  & IN  METHOD  = ’ CORRECT'  , (Description  of  2nd 
gas  by  2nd  path  path)  4 END. 


Wrong  adjustment  of 
hot  gas  by  3rd  path 


4 IN  METHOD  = ' WRONG*  , (Description  of  3rd 
path)  4 END 


It  is  emphasized  that,  if  METHOD  is  not  specified,  the  current  namelist  is 
taken  as  a hot  gas,  as  opposed  to  an  intervening  atmospheric  path,  and  subsequent 
namelistswith  METHOD  specified  will  be  treated  as  atmospheric  paths  applied  to  it. 
Note  that  the  first  namelist  executed  should  describe  a hot  gas  or  other  source  ; 
specification  of  METHOD  in  the  first  namelist  will  lead  to  erroneous  (and  probably 
disastrous)  results . 

Ill)  The  third  change  made  in  the  program  is  that  the  determination  of  the  expo- 


nent in  the  overlap  factor  (/?/ B )^  , described  on  page  4 of  the  28  January  memo, 


has  been  changed  so  that  certain  extreme  paths  are  treated  correctly.  The  form 
that  rj  now  takes  is  quite  complex  and  is  described  in  a separate  memo  to  be  written 
shortly.  The  following  paragraph  explains  the  options  available  and  should  replace 
the  paragraph  beginning  in  the  middle  of  page  4 of  the  memo  of  January  28,  1970. 

The  exponent  17  of  the  overlap  factor  (/?/ pJP  , can  be  treated  in  two  ways 

1)  set  equal  to  namelist  input,  EXPO 

2)  computed  according  to  the  approximation  described  in  a later 
memo  on  this  subject. 

Also,  the  description  of  the  use  of  EXPO  on  page  7 of  the  28  January  memo  should  be 
changed  to  read 


EXPO 


{ 0.  * EXPO  * 1.,  EXPO 

> *-  °r  <°-  1 computed  from  approximation  described  in 
^ later  memo . 
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THE  UNIVERSITY  OF  MICHIGAN 
THE  INFRARED  AND  OPTICS  LABORATORY 


MEMO  TO:  3211  File 

FROM:  G.  Lindquist^^/^ 


22  July  1970 


SUBJECT:  How  to  use  the  program  PLUMES 

PLUMES  is  a program  designed  to  compute  the  radiant  intensity  from  a 
cloud  of  hot  gas . This  is  done  by  dividing  the  projected  area  of  the  cloud  into 
elements,  and  computing  the  temperature  and  pressure  distribution  along  a line 
of  sight  through  the  cloud  at  the  center  of  each  element.  A band  model  calcula- 
tion for  the  temperature  and  pressure  distribution  along  a particular  line  of 
sight  done  in  a manner  identical  to  that  performed  by  the  program  FS  (Ref.  1)  - 
yields  the  radiance  along  that  line  of  sight.  The  summation  of  the  product  of 
such  radiances  times  the  projected  area  of  the  element  associated  with  that  line 
of  sight  yields  the  radiant  intensity  of  die  total  gas  cloud  as  viewed  at  a distance. 

In  the  current  formulation  the  program  was  designed  to  handle  rocket  and  jet  en- 
gme  exhaust  plumes  so  that  its  use  is  limited  to  axisymmetric  sources . Sufficient 
generality  is  included  so  the  gas  cloud  can  be  viewed  from  any  aspect.  Also  an 
atmospheric  path  can  optionally  be  inserted  between  the  gas  cloud  and  the  observer. 

Two  absorbing  gases  and  a broadening  gas  can  be  included.  Normally  these  are 
H20,  C02andN2. 


i 


The  number  of  input  parameters  is  large  and  they  can  best  be  described 
in  groups  corresponding  to  their  function  in  the  program.  All  the  parameters 

escribed  here  are  read  in  by  means  of  a namelist  statement  from  input  unit  1 . 

The  name  of  the  namelist  is  IN. 

|$  * 

* % i 
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Plume  Specification 

The  gas  cloud  is  assumed  to  be  axisymmetric  so  that  it  is  entirely  specified 
by  specifying  the  conditions  on  one  side  of  a plane  containing  the  plume  axis.  The 
plume  is  so  specified  by  specifying  the  temperature  and  total  pressure  distributions 
in  such  a plane  using  polar  coordinates  with  origin  at  the  center  of  the  nozzle  exit 
plane  (all  of  the  plume  is  taken  to  lie  downstream  of  the  nozzle  exit  plane). 

This  specification  can  be  handled  in  two  ways,  controlled  by  the  variable 
SETUP3.  If  SETUP3  is  negative  the  plume  parameters  are  read  from  input  unit  3 
according  to  a format  described  in  Appendix  A.  In  this  mode,  the  composition  of 
the  plume  is  fixed  in  the  following  way.  The  inputs  from  unit  3 are  only  temperature 
and  total  pressure.  The  partial  pressures  of  the  three  gases  used  are  then  deter- 
mined by  the  parameters  PARP1,  PARP2,  and  PARP3  which  are  die  mole  fractions  of 
gases  1,  2,  and  3 respectively.  Options  are  also  available  to  allow  the  input  param- 
eters read  in  on  unit  3 to  be  dimensionless  temperatures,  pressures,  and  distances. 
These  options  are  contained  in  the  variables  PZERO,  TZZRO,  and  RRATIO.  The 
values  read  in  on  unit  3 for  temperature  are  multiplied  by  TZZRO,  for  pressure  by 
PZERO  and  for  distance  by  RRATIO.  If  any  of  these  options  are  not  to  be  used,  the 
particular  parameter  should  be  set  equal  to  1 . 

If  SETUP3  is  set  positive  it  is  assumed  that  the  plume  specifications  have 
already  been  read  from  unit  3 in  a previous  calculation. 

If  SETUP3  is  zero,  the  plume  parameters  are  calculated  using  an  idealized 
gas  dynamic  plume  model  (see  Ref.  2)  characterized  by  an  initial  region,  in  which 
mixing  occurs  adjacent  to  a potential  core,  and  a main  region  in  which  mixing 
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with  the  ambient  determines  the  physical  parameters . Figure  1 shows  a dia- 
gram of  such  a plume . The  boundary  of  this  plume  model  also  determines  the 
extent  of  the  area  over  which  the  plume  integration  takes  place.  Thus,  it  is 
necessary  to  include  the  quantities  which  describe  this  boundary  even  though 


SETXJP3  is  negative, 

i.e.  plume  pressures  and  temperatures  are  to  be  read  in 

rather  than  calculated . 

The  parameters  which  describe  the  gas  dynamic  plume  model  are  as 
follows: 

RO 

the  radius  of  the  nozzle  exit  (cm) 

ZC 

the  length  of  the  mixing  region  of  the  plume  (cm) 

PANGD 

the  angle  between  the  plume  axis  and  the  outer  boundary 
of  the  plume  in  the  mixing  region  (degrees) 

ZMAX 

the  total  length  of  the  plume  (cm) 

zee 

the  length  of  the  potential  core  of  the  plume  (cm) 

If  SETUP3  is  zero, 

parameters  are  needed  which  describe  the  physical  properties  of 

the  gas  dynamic  model.  These  are  described  as  follows: 

TPL1 

temperature  in  the  potential  core 

TPL2 

ambient  temperature  (temperature  of  gas  surrounding  plume) 

Pll 

pressure  of  gas  1 is  the  potential  core 

P12 

ambient  pressure  of  gas  1 

P21 

pressure  of  gas  2 in  the  potential  core 

P22 

ambient  pressure  of  gas  2 

P31 

pressure  of  gas  3 in  the  potential  core 

P32 

ambient  pressure  of  gas  3 

7* 
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To  repeat,  RO,  ZC,  PANGD  and  ZMAX  must  be  included  as  input  even  if  SETUP3 
is  not  zero  to  define  the  region  over  which  integration  takes  place.  The  remainder 
of  the  parameters  describing  the  gas  dynamic  plume  model  are  needed  only  if  SETLJP3 
= 0. 

Atmospheric  Path  Specification 

In  order  to  include  an  atmospheric  path  in  front  of  the  plume  it  is  only  neces- 
sary to  specify  die  altitude  of  the  plume  and  of  the  observer  in  cm  (ALTPLM  and 
ALTOBS  respectively),  the  range  between  the  plume  and  the  observer  in  cm  (RANGE) 
and  the  number  of  elements  into  which  the  atmospheric  path  is  to  be  broken(NA).  The 
proper  temperatures  and  pressures  are  obtained  from  a model  atmosphere  contained 
in  the  program.  If  no  atmospheric  path  is  to  be  inserted  set  RANGE  = 0. 

Frequency  Range 

It  is  possible  to  includetwo  groups  of  frequencies  in  any  one  calculation.  The 
first  group  is  specified  by  the  variables  NFREQ,  NUFRST,  and  NUINC.  NUFRST 
describes  the  first  frequency  to  be  computed,  NUINC  describes  the  spacing  between 
frequencies  and  NFREQ  describes  the  number  of  frequencies  in  the  first  group  of  fre- 
quencies . The  second  group  of  frequencies  is  described  by  the  variables  NFREQ1, 
NUFRS1,  and  NUINC  1 which  correspond  exactly  in  usage  to  NFREQ,  NUFRST  and 
NUINC.  Thus,  for  example  in  one  calculation  the  frequency  ranges  from  2100  to 
3000  and  from  3500  to  4000  cm’1  can  both  be  covered.  If  only  one  group  of  frequencies 
is  to  be  covered,  NFREQ1,  NUFRS1,  and  NUINC1  can  be  ignored.  Alternatively, 
NUFRST  can  be  set  negative  and  the  program  then  expects  to  find  NFREQ  values  of 
frequency  read  in  as  the  array  NU. 
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Bandmodel  Parameters 

TTie  bandmodel  parameters  for  the  absorbing  gases  are  expected  to  be  in 
external  files  in  the  expected  order.  The  order  is  the  same  as  was  given  for  use 
in  the  program  FS  (see  Ref.  1).  The  band  model  parameters  for  the  first  absorb- 
ing gas  (usually  H20)  are  read  from  input  unit  4,  and  the  band  model  parameters 
for  the  2nd  absorbing  gas  'usually  C02)  are  read  from  unit  5.  The  array  GAS 

contains  the  names  of  the  gases  to  be  considered,  usually  'H20',  'C02'  , and  'N2' 
in  that  order . 

Fineness  of  the  Plume  Integration 

Integration  takes  place  In  three  dimensions  in  this  program.  The  plume 
is  projected  into  a plane  which  is  perpendicular  to  the  line  of  sight  from  the  plume 
to  the  observer  and  passes  through  the  center  of  the  nozzle  exit.  Two  dimensions 
of  the  integration  occur  as  integration  over  the  area  of  the  plume  as  projected  into 
this  plane.  The  origin  of  the  coordinate  system  used  in  this  plane  is  at  die  center 
of  the  nozzle  exit.  Polar  coordinates  are  used.  The  third  dimension  is  perpendicu- 
lar to  the  above  mentioned  plane  and  occurs  as  integration  along  each  of  the  lines 
of  sight  through  the  plume . To  obain  a line  of  sight  radiance  as  the  result  of  this 

integration,  this  integration  is  carried  out  over  the  transmitance  along  the  line  of 
sight. 

The  fineness  of  the  transmitance  integration  is  controlled  by  a variable 
NMAX.  The  maximum  and  minimum  expected  temperatures  are  read  in  as  TMAX 
and  TMIN  and  this  range  is  divided  into  NMAX  temperature  increments . In  calcu- 
lating through  a line  of  sight , each  time  die  boundary  of  one  of  these  temperature 
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V increments  is  crossed,  a transmittance  element  is  established.  Thus,  a line  of 


sight  which  traverses  the  full  range  of  TMAX  to  TMIN  will  be  divided  into  NMAX 
elements . Lines  of  sight  with  alternately  increasing  and  decreasing  temperatures 
will  be  divided  into  more  elements,  and  lines  of  sight  having  only  small  tempera- 
ture variations  will  be  divided  into  correspondingly  fewer  elements.  The  tempera- 
ture is  divided  into  increments  according  to  the  inverse  2nd  power  of  the  temperature 
so  that  higher  temperatures  yield  smaller  increments . This  procedure  is  useful 
without  modification  for  plumes  having  only  a limited  temperature  range . If  the 
temperature  range  is  large,  or  if  the  temperatures  are  clustered  around  two  or 
more  values,  such  a breakup  is  not  very  advantageous . For  such  occurrences, 
provisions  have  been  made  to  allow  the  temperatures  defining  the  temperature  in- 
crements to  be  read  in  directly.  This  is  done  by  setting  ICHOOZ=i  and  reading  in 
NMAX+3  values  for  the  temperature  increment  boundaries  using  the  array  TLIM, 
e.g.  by  including  in  the  namelist  IN  a line  such  as  ICHOOZ=l,  TLIM  = t.,  t , --- 

tNMAX+3’  where  ti  tNMAX+3  416  t*le  temPerature  increment  boundaries.  If  the 
values  of  TLIM  are  read  in  TMAX  and  TMIN  can  be  ignored. 

The  fineness  of  the  integration  over  the  radial  dimension  in  the  plane  of  pro  - 
jection  is  also  controlled  by  NMAX  (the  controlling  of  integration  fineness  in  two 
dimensions  by  one  variable  resulted  in  some  loss  of  flexibility  but  was  done  to  reduce 
die  number  of  input  variables  that  needed  to  be  considered) . The  value  of  the  radial 
increment  is  chosen  so  that,  with  a blob  of  ga3  which  has  equal  extent  in  all  directions, 
the  number  of  radial  increments  will  approximately  equal  the  average  number  of  line 
of  sight  transmittance  increments.  In  the  edge  regions  of  the  plume  the  radial  incre- 
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ment  size  is  allowed  to  increase  with  decreasing  average  temperature.  The  reason 
for  this  is  to  reduce  the  computing  time  by  computing  fewer  radiance  values  for 
regions  where  the  radiance  will  be  very  samll  compared  to  the  peak  radiance  values. 

The  fineness  of  the  integration  over  the  angular  dimension  in  the  projected 
plane  is  not  variable.  This  angular  dimension  is  called  PHI  in  the  program . This 
angle  for  a given  line  of  sight  lies  in  the  plane  of  projection.  It  is  the  angle  between 
the  projection  of  the  plume  axis  into  the  above  mentioned  plane  and  the  line  from  the 
center  of  the  nozzle  exit  to  the  intersection  of  the  given  line  of  sight  with  the  same 
Plane  (all  lines  of  sight  through  the  plume  are  taken  to  be  parallel  to  the  line  of  sight 
from  the  observer  to  the  plume  so  that  all  lines  of  sight  are  perpendicular  to  the 

above  mentioned  plane).  The  values  of  PHI  currently  used  in  the  integration  are:  0°, 

6 , 18  , 36°,  60°,  90°,  120°,  144°,  162°,  174°,  180°.  Integration  over  180°  < PHI 
< 360  is  accounted  for  by  multiplying  the  result  for  0 < PHI  < 180°  by  two . 

Modifications  to  the  above  integration  procedure  can  be  made  by  simple  modi- 
fications to  the  program  itself. 

Other  Features 

Three  auxiliary  temperatures  are  available  for  use.  If  TBaCK  is  set  non-zero, 
a black  background  at  temperature  TBACK  will  be  placed  behind  the  plume.  If  TBACKC 
is  set  non -zero,  radiance  values  equal  to  a black  body  at  TBACKC  will  be  subtracted 
from  all  radiance  values  as  a background  compensation.  If  TOOZ  is  non-zero,  in  case 
where  foe  line  of  sight  looks  upthenozzlea  black  background  at  TNOZ  will  be  used  to 
represent  the  inside  of  die  nozzle.  In  addition  if  OBSCUR  is  set  non-zero,  an  opaque 
disc  of  radius  OBSCUR  in  cm  will  be  placed  a,  the  nozzle  exit  and  any  lines  of  sight 
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which  intersect  this  disc  will  be  ignored. 


Following  i^a  glossary  of  all  variables  which  can  be  set  via  namelist  IN. 
Not  all  are  mentioned  in  the  text  since  they  are  merely  control  parameters . Also 
a group  of  parameters  controls  the  scales  used  in  making  plotter  cards  for  the 
Benson -Lehner  plotter.  The  details  of  the  transmittance  calculation  are  exactly 
as  described  in  Ref.  1. 


The  program  can  be  run  by  the  use  of  the  following  statement  (assuming 
it  is  run  from  SX90)  $RUN  PO+SX27:SO  + ST82:SUBPACK+NTAU50BJ;  l=<data> 
3=<plume  data>  4=BKH20  5=BKC02  9=  <plotter  cards> 

The  execution  time  of  this  program  is  long  and  it  occupies  a great  deal  of  computer 
memory  so  that  it  is  economically  feasible  to  run  it  only  in  batch  mode.  In  the  above 
statement  words  enclosed  in  M<  >"  e.g.  <data>  indicate  filenames  where  that 
particular  information  either  is  scored  on  input,  or  to  be  stored  on  output.  Printed 
output  comes  from  unit^/and  unit  8 . 
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GLOSSARY  OF  INPUT  PARAMETERS 


Data  is  read  in  by  means  of  namelist  IN,  the  elements  of  which  are: 


SETUP3 


SETUP1 


SETUP2 


NUFRST 

NFREQ 

NUINC 

NU 

NUFRS1 
NFREQ1 
NUINC  1 
NA 

ALTOBS 


> -0,  plume  conditions  are  not  read  in 

= 0.,  the  plume  conditions  are  claculated  from  certain  required 
input  data  (TPL1,  TPL2,  Pll,  P12,  P21,  P22,  P31,  P32) 

< -1.,  plume  conditions  read  in  by  means  of  3=FDNAME  in  $RUN 

instructions . FDNAME  is  the  file  in  which  the  information 
is  stored. 

> 1.,  k,  j8  for  absorbing  gas  1 not  read  in. 

< 1.,  k,  jS  for  absorbing  gas  1 read  in  by  means 

of  4=FDNAME  in  $RUN  instructions. 

FDNAME  is  the  file  where  the  data  is  stored. 

> 1 . , k , p for  absorbing  gas  2 not  read  in . 

< 1 . , k,  /3  for  absorbing  gas  2 read  in  by  means 

of  5=FDNAME  in  $RUN  instructions. 

FDNAME  is  the  file  where  data  is  stored. 

first  wavenumber  in  first  group  of  wavenumbers  for  which  calculations 
are  to  be  made  (real  number) . 

=0. , wavenumbers  are  not  to  be  calculated  but  to  be  read  in  via  NU 

number  of  wavenumbers  in  first  group  of  wavenumbers  for  which  the 
calculations  are  to  be  made . 

(integer) 

wavenumber  increment  (in  wavenumbers)  in  first  group  of  wavenumbers 
(If  NUFRST=0.,  NUINC  can  be  ignored) 

wavenumbers  (real  numbers)  (If  NUFRST=0.,  the  wavenumbers  are 
read  in;  if  not,  NU  can  be  ignored) 

first  wavenumber  in  2nd  group  of  wavenumbers . 

number  of  wavenumbers  in  2nd  group  of  wavenumbers  (integer) 

wavenumber  increment  in  2nd  group  of  wavenumbers . 

number  of  increments  for  the  atmospheric  path  (integer  s 20) 

altitude  of  the  observer  or  observing  instrument  (real  number,  in  cm) 


These  variables  can 
almost  always  be  ignored 
since  the  program  sets 
them  automatically. 
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ALTPLM 

RANGE 

ALPHA 

GAS 

IATMO 

RO 

PANGD 

ZC 

ZCC 

OBSCUR 

ZMAX 

TMAX 

TMIN 

ASPDEG 

NMAX 

TBACK 

TBACKC 

TN02 


altitude  of  the  exhaust  plume  (real  number,  in  cm) 

distance  from  observer  to  plume  (real  number,  in  cm) 

broadening  coefficients  for  absorbing  gases  1 and  2.  ( .2  for  H20-N2, 
and  .769  for  C02-N2)  Can  be  ignored  of  the  gases  are  'H20'  'C02' 
and  'N2'. 

names  of  the  two  absorbing  gases  and  the  broadening  gas . Can  be  ignored 
if  the  gases  are  to  be  * H20'  , 'C02'  and  ' N2'  . 

> 0,  atmospheric  path  information  is  printed. 

< 0,  atmospheric  path  information  is  not  printed. 

radius  of  plume  at  nozzle  exit  (real  number,  in  cm) 

plume  angle  (real  number,  in  degrees) 

length  of  plume  after  which  plume  diameter  remains  constant . 

(real  number, in  cm) 

length  of  potential  core  of  plume  (real  number,  cm). 

obscuration  of  line  of  sight  at  s=0.  (Radius  of  circular  obscuration 
at  nozzle  exit  centered  about  plume  axis)  (real  number,  in  cm). 

maximum  length  of  plume 

maximum  temperature  of  plume  (°K) 

minimum  temperature  of  plume  (°K) 

aspect  angle  in  degrees  (real  number) 

approximate  number  of  increments  into  which  two  of  the  three  dimen- 
sions is  divided.  The  total  number  of  cubes  into  which  the  plume 
will  be  divided  a:(NMAX)  -11 

temperature  of  background  behind  plume  (°K) 

temperature  of  background  radiance  which  is  to  be  subtracted  from 
radiance  values  ( K) 

temperature  of  inside  of  nozzle  (used  if  line  of  sight  looks  up  the 
nozzle)  K.  r 
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TODAY  date  program  is  run  (format  real  *A8) 

ITYPE  1,  overlapping  lines 

2,  nonoverlapping  lines 


EXPO 

INPRNT 

IPLOT 

IDE 

XMAX 

XMIN 

ISCALE 

RMAX 

LCODE 

ITAU 

IL 


0 ^EXPO  s \ ^etermine^  ^ sPecial  interpolation  procedure  (Ref.  3) 

> 0,  path  information  printed 

< 0,  path  information  not  printed 

< 0,  results  not  punched  on  cards  for  Benson  Lehner  plotter 

> 0,  results  punched  on  cards  for  IPLOT  plots . 


label  for  punched  data  (format  3A4) 


array  giving  maximum  abscissa  value  for  each  of  the  IPLOT  plots 
(wavenumbers)  (real  number)  (should  have  IPLOT  values) 

array  giving  minimum  abscissa  value  for  each  of  the  IPLOT  plots 
(wavenumbers)  (real  number)  (should  have  IPLOT  values). 


These  varia- 
bles control 
the  output  fo: 
the  Benson 
Lehner  plottt 


< 0,  automatic  scaling  off,  maximum  ordinate  value  for  plots 
taken  from  RMAX 
> 0,  automatic  scaling  on. 

array  giving  maximum  ordinate  value  (radiant  intensity)  (real 
numbers)  for  each  of  the  IPLOT  plots  (should  have  IPLOT 
values  if  ISCALE  < 0). 

6,  points 
0,  lines 


these  variable 
control  the  ou 
> put  for  the 
Benson  Lehnei 
plotter. 


<0, 

>0, 

<0, 

>0, 


transmittances  are  not  printed 

transmittances  are  printed  for  ITAU  wavenumbers. 

line  of  sight  radiances  are  not  printed 

line  of  sight  radiances  are  printed  for  IL  wavenumbers 


If  data  for  plume  is  to  be  calculated  (SE7UP3=0.  in  namelist),  the  following 
elements  must  be  in  namelist;  otherwise  they  can  be  ignored. 

TPL1  temperature  in  core  (real  number,  °K) 

TPL2  temperature  at  edge  (real  number,  °K) 
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PH  pressure  of  GAS(1)  in  core  (real  number,  atm) 

P!2  pressure  of  GAS(1)  at  edge  (real  number,  atm) 

P21  pressure  of  GAS(  2)  in  core  (real  number,  atm) 

P22  pressure  of  GAS(2)  at  edge  (real  number,  atm) 

P31  pressure  of  GAS(3)  in  core  (real  number,  atm) 

P32  pressure  of  GAS(3)  at  edge  (rerJ  number,  atm) 

Some  of  the  namelist  elements  are  preset  by  means  of  DATA  statements: 


SETUP3=  -10. 


EXPO  = -9, 


SETUP1=-10. 


SETUPS  -10. 


NUFRST  = 0 . 

NA  = 10 
ALPHA(1)  = .2 
ALPHA(2)  = .769 


IATMO  = -1 


IPLOT  = -1 


IDE  = 'RAD  VS  NU' 


LCODE = 6 


ITAU  = -1 
IPRNT  = -1 


IL  = -1 


I TYPE  = 1 


TNOZ=  900.  (°K) 

TBACK  = 0. 

TBACKC  = 0 . 

GAS  = 'H20'  , 'C02' , 'N2' 
ISCALE  = 1 


J n>i.  •imi  ima 
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Extracted  From: 


IR  Radiation  From  Rocket  Exhaust  Plumes  in  Space1*,  Final  Report 
17  July  67  - 7 July  70.  F.  Simons  and  G.  Lindquist,  Oct  70,  Report 
3211-1-F  contract  F04701-70C-0004. 
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METHODS  FOR  THE  CALCULATION  OF  INFRARED  RADIATION  FROM  EXHAUST  PLUMES 

<U)  Aside  from  condensed  particle  emission,  the  infrared  radiation  from  an  exhaust  plume 
consists  of  the  spectral-line  emission  in  the  vibration -rotation  bands  of  the  optically  active 
products  of  combustion.  For  most  operational  rocket  engines,  excluding  those  utilizing  fluorine 
or  it.  compounds  as  oxidizers,  carbon  dioxide  and  water  vapor  are  the  dominant  radiators.  For 


these  well-known  molecules,  modern  quantum -mechanical  methods  can  easily  provide  the  line 
positions  corresponding  to  the  differences  in  the  vibrational  and  rotational  energy  levels  asso- 
ciated with  emission  or  absorption  of  photons.  However,  at  higher  temperatures,  many  vibra- 
tional and  rotational  levels  become  populated,  and  the  number  of  spectral  lines  contributing  to 
the  radiation  increases  enormously.  There  is  not  presently,  nor  will  there  be  in  the  foreseeable 
future,  sufficient  data  on  the  line  strengths  and  widths  to  perform  exact  calculations  of  radiative 
transfer  for  such  hot  gases.  Moreover,  even  if  these  data  were  available,  their  application  In  a 
line  by  line  calculation  of  the  spectral -radiance  distribution  and  radiant  intensity  in  a complex 
flow  field  would  be  quite  impractical,  even  with  modern  high-speed  computers. 

(U)  Consequently,  approximate  methods  must  be  Invoked  to  handle  the  treatment  of  rrdiation 
from  highly  nonuniform  and  nonlsothermal  gases  such  as  exhaust  jets.  The  most  sophisticated 
of  these  approximate  methods  for  treating  infrared  radiative  transfer  Is  the  use  of  a molecular 
band  model  to  obviate  the  requirement  of  considering  the  variation  in  absorption  and  emission 
in  two  dimensions:  spatial —along  the  nonlsothermal  path  through  the  radiating  region;  and 
spectral -over  the  contour  of  each  spectral  line.  In  essence,  a molecular  band  model  is  a 
simple  mathematical  representation  for  the  location  and  distribution  in  strengths  of  the  rotation- 
al lines:  the  most  generally  useful  model  is  one  which  assumes  a random  location  of  lines  (an 
assumption  which  is  quite  realistic  for  polyatomic  molecules  in  which  many  vibrational  and 
rotational  energy  levels  become  populated  in  several  vibrational  degrees  of  freedom).  By  this 
means,  the  consideration  of  the  variation  in  emission  and  absorption  over  the  contours  of  in- 
dividual spectral  lines  can  be  obviated,  and  radiative  transfer  expressions  developed  which  re- 
quire only  consideration  of  variations  along  the  optical  path.  Thus,  spectral  radiances  can  be 
calculated  over  the  projected  area  of  a radiating  body  of  gas  (e.g.,  an  exhaust  plume)  and  inte- 
grated to  yield  spectral  radiant  Intensities. 

(U)  Band  models  are  generally  formulated  in  terms  of  two  parameters  which  are  combina- 
tions of  the  average  values  of  three  basic  variables:  the  line  strengths,  widths,  and  separations. 
These  parameters  are  functions  of  both  wavelength  and  temperature  for  band  models  applied  to 
gases  which  emit  as  well  as  absorb,  and  are  determined  experimentally  in  the  laboratory  from 
measurements  of  absorption  In  isothermal  samples  of  gas.  Their  use  in  calculations  of  radiant 
emission  for  nonlsothermal  paths  then  Involves  the  above  radiative -transfer  relations. 

(U)  Molecular  band  models  were  originally  developed  some  years  ago  for  the  treatment  of 
atmospheric  absorption  in  meteorological  applications;  a review  of  the  state  of  the  art  in  this 
area  has  been  published  as  a Willow  Run  report  [10).  The  extension  of  band  models  to  handle 
the  case  of  inhomogeneous  hot  gases  has  been  the  subject  of  recent  research  at  The  University 
of  Michigan  [11-18] , UCLA  [ 19] , General  Dynamics/Conva'r  [20,22],  and  Warner  and  Swasey 
[21]  • For  the  most  part  the  motivation  for  these  studies  was  the  calculation  of  radiation  from 
rocket  exhaust  plumes:  this  subject  was  reviewed  at  a recent  NASA -sponsored  symposium  [17].' 
The  inhomogeneous  gas  band-model  expressions  developed  at  The  University  of  Michigan  repre- 
sent a further  advance  beyond  the  nonlsothermal  band -model  formulations  developed  elsewhere 
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and  based  on  the  Curtis -Godson  approximation.*  The  University  of  Michigan  band  models  in- 
volve application  of  the  "nearly -weak"  line  and  "nearly -strong"  line  approximations  [12,  14] , 
each  of  which  reduces  to  the  Curtis -Godson  approximation  and,  thus  each  provides  an  identical 
value  when  applied  to  the  calculation  of  transmission  through  an  inhomogeneous  gas.  However, 
when  applied  to  the  calculation  of  emission,  these  two  approximations  provide  significantly 
different  results,  especially  for  the  case  of  optical  paths  exhibiting  a large  variation  in  tempera- 
ture. Accordingly,  each  has  an  appropriate  range  of  optical  depths,  and,  for  the  intermediate 
region,  an  interpolation  between  the  two  can  be  made. 


(U)  The  resultant  formulation  for  this  band  model  in  terms  of  the  average  spectral  radiance 
(specific  intensity)!  , the  Planck  function  1*(T)  of  the  local  temperature,  and  average  spectral 


absorptance  a(u)  in  the  frequency  interval  is  [12,  18] 


ft(Xr  ) „ 

V^eJ  ^“P  <-<y>df 


a(y)  = 1 - exp 


h«v] 


The  function  f is  the  Laienburg-Reiche  function  (or  a simple  algebraic  approximation  thereof) 
of  the  dimensionless  optical  depth  x,  in  accordance  with  the  following  definitions: 


f ■ f(x)  = xe'x[j0(i£)  - iJjfbdJ  = x/(l  ♦ nx/Z)l/2 

x = p'1  fV/fYdX' 

6 JO  VV 


The  quantity  k = k(i>,  T)  is  the  first  band-model  parameter:  the  average  absorption  coefficient 
which  can  be  identified  with  the  ratio  of  average  line  strength  to  spacing.  The  quantity  p = p(v,  T) 
is  the  second  band-model  parameter:  the  line  overlap  factor  which  is  proportional  to  the  ratio 


of  line  width  to  spacing.  The  effective  value,  0e,  Is  the  weighted  average  along  the  entire  path 


of  length  L: 


f^kdX 


f^kdX 


where  X is  the  optical  depth  in  centimeter -atmospheres  normalized  to  standard  conditions.  The 
exponent  rj  provides  for  the  interpolation  between  the  nearly  weak  and  nearly  strong  line  approxi- 
mations by  means  of  the  empirically  derived  specif ications:  x L < 1,  77  = 0;  1 s xL  5 10 , rj  = 


*(U)  The  Curtis -Godson  approximation,  in  essence,  provides  the  means  for  representing 
an  Inhomogeneous  path  by  an  equivalent  homogeneous  one:  it  was  originally  developed  to  facili- 
tate calculations  of  absorption  along  nonuniform  paths  in  the  atmosphere. 
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[(xL  - !)/•) ; xL  > 10,  i?  = 1.  The  validity  of  this  band  model  for  calculation*  of  radiation  from 
inhomogeneous  bodies  of  hot  gases  is  indicated  in  Figs.  2 and  3 ; these  figures  show  predicted  and 
observed  spectra  for  samples  of  H20  and  C02  at  a pressure  of  1 atm  in  a 60  cm  path  with  a 

symmetrical  triangular  profile  in  temperature  with  values  of  1250°K  at  the  center  and  400°K  at 
the  ends  [18]. 

(U)  The  preceding  examples  Illustrate  the  capability  of  this  band  model  to  provide  reliable 
values  of  infrared  radiances,  both  in  magnitude  and  spectral  distribution,  for  moderate  optical 
depths  of  water  vapor  and  carbon  dioxide  in  which  temperatures  vary  by  a factor  of  3 or  so.  A 
more  severe  test  of  the  model  is  an  application  to  a region  In  which  the  temperature  varies  by 
an  order  of  magnitude.  To  test  the  model  for  this  case,  a number  of  measurements  was  re- 
cently made  at  the  Willow  Run  Laboratory  (under  a NASA -supported  study  of  Infrared  spectros- 
copy as  a diagnostic  tool  [541  of  the  infrared  emission  from  2.5 -cm  methane -oxygen  flame. 

The  results  are  shown  in  Fig.  4.  The  upper  spectrum  is  that  of  the  flame  viewed  directly-  the 
smoother  line  represents  the  results  of  the  band-model  calculation  based  on  the  flame  model, 
shown  in  Fig.  5,  which  takes  into  account  variations  in  HzO  and  C02  concentrations  due  to  both 
recombination  and  dilution  in  the  mixing  region.  The  temperature  in  the  flame  varied  from 
tne  adiabatic  flame  value  of  3000°K  at  the  center  to  the  ambient  value  of  300°K  at  the 
boundary;  this  close  agreement  speaks  well  for  the  band  model.  However,  a limitation  of  this 
model  appeared  in  applications  involving  a small,  hot  source  Imbedded  in  a large,  cool,  absorb- 
ing region.  The  lower  spectrum  In  Fig.  4 represents  such  a case;  the  2.5-cm  flame  viewed 
through  a 10-m  path  of  ambient  air.  Here  the  band  model,  In  the  form  specified  above,  exhibited 
poor  agreement  and  anomalous  behavior,  and, hence,  required  some  modification.  The  modified 
band  model  then  yielded  predictions  in  close  agreement  with  the  observations,  as  indicated  in  " 
Fig.  4.  The  modification  consists  of  an  Improved  method  of  interpolation  between  the  nearly 
weak  and  nearly  strong  line  approximations,  i.e.,  the  evaluation  of  V,  which  was  derived  from  a 
comparison  with  calculations  using  the  exact  equation  of  transfer.  This  subject  will  be  dis- 
cussed in  detail  in  another  publication.  For  purposes  of  the  present  study,  the  differences  in 
the  predictions  of  the  band  model  between  the  original  and  modified  form  are  negligible. 
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THE  PLUME-RADIATION  COMPUTER  CODE 

(U)  A computer  code  has  been  compiled  for  the  calculation  of  spectral-radiance  distribu- 
tions and  spectral-radiant  intensities  of  exhaust  plumes  viewed  from  long  ranges  at  various 
aspects.  The  main  part  of  the  code  consists  of  a program  designated  PLUMES,  in  which  plume 
properties  (i.e.,  temperatures  and  species  partial  pressures)  are  put  in  as  functions  of  radial 
and  axial  distances;  therefore,  at  the  moment,  the  code  is  limited  to  axlsymmetric  sources. 

T*e  aspect  angle  is  specified,  and,  through  a transformation  of  coordinates,  the  plume  proper- 
ties are  determined  along  the  various  parallel  lines  of  sight  through  the  plume  over  the  entire 
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FIGURE  6 DIAGRAM  OF  THE  TWO  COORDINATE  SYSTEMS  USED  IN  INTEGRATING  OVER 

THE  PLUME 

projected  area  in  the  line  of  sight,  as  indicated  in  Fig.  6.  Provisions  are  made  for  varying  the 
Incremental  sizes,  both  of  the  elements  in  the  plume  projection  and  the  linear  segments  of  the 
optical  path  through  the  gas,  with  the  finer  resolution  being  near  the  nozzle  exit  where  the  steep- 
est gradients  occur.  These  values  then  become  input  to  the  band- model  subroutine,  which  has 
been  designated  XTAU. 

(U)  The  band  model,  as  specified  in  Eqs.  (1)  through  (5),  has  been  programmed  by  means 
of  the  trapezoidal  approximation.  The  band-model  parameters,  extracted  from  the  extensive 
tabulations  of  the  GD/Convalr  group  [20] , are  contained  in  two-dimensional  files  as  functions 
of  temperature  and  frequency.  An  initial,  final,  and  Increment  of  frequency  are  specified. 

(U)  The  calculation  then  proceeds  as  follows.  The  PLUMES  program  starts  at  the  first  in- 
crement of  area  on  the  projection  of  the  plume  and  specifies  the  conditions  along  the  path  through 
the  gas.  The  XTAU  subroutine  then  starts  at  the  first  frequency,  calculates  the  net  spectral 
radiance  at  the  boundary,  proceeds  through  successive  frequencies  until  the  last,  and  returns 
the  values  of  spectral  radiance  vs.  frequency  to  PLUMES  for  storage.  PLUMES  then  submits 
the  path  conditions  for  the  second  Incremental  area,  and  the  process  is  repeated  until  the  entire 
area  of  the  plume  projection  la  covered. 
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The  spectral  radiances  at  the  specified  frequencies  are  then 
summed  to  yield  spectral- radial  intensities.  The  program  then  proceeds  to  the  next  aspect 
angle  and  repeats  the  process.  In  addition  to  the  primary  output  of  spectral-radiant  intensities 
as  functions  of  frequency  and  aspect  angle,  spectral-radiance  maps  at  specified  frequencies 
and  angles  can  be  called  for. 

(U)  A few  details  of  the  spectral-radiance  calculations  are  as  follows.  The  nonisothermal 
path  through  the  plume  is  described  by  a number  of  isothermal-path  elements  of  varying  length. 

At  a particular  frequency,  the  radiance  is  computed  by  evaluating  the  radlative-tranafer  equa- 
tion in  the  following  form 

I ./Vdi  (6) 

" Jl  * 

which  is  simply  an  alternative  form  of  Eq.  (1).  Thus,  it  is  necessary  to  determine  both  I*  (the 
blackbody  function)  and  r as  a function  of  position  along  the  path.  The  routine  XTAU  computes 
the  transmittance  from  the  beginning  of  the  path  to  each  of  the  boundaries  between  path  elements. 
These  values,  in  conjunction  with  the  blackbody  function,  are  used  in  evaluating  Eq.  (6)  by  the 
trapezoidal  rule. 

(U)  Upon  the  choice  of  a frequency  and  a set  of  elemental  path  lengths,  temperatures,  and 
pressures,  the  routine  XTAU  first  computes  0e  by  evaluating  Eq.  (5).  Then,  x,  given  by  Eq.  (4), 
is  computed  from  the  beginning  of  the  path  to  each  of  boundaries  between  path  elements.  The 
revised  Interpolation  technique  mentioned  earlier  is  used  to  determine  jj.  In  this  revised  pro- 
cedure, i)  is  Itself  a function  of  x (as  well  as  of  £/0e).  Thus,  the  exact  determination  of  x in- 
volves a trial  and  error  procedure  in  which  a value  of  h is  estimated,  a value  of  x then  calcu- 
lated, and,  t),  in  turn,  recalculated  until  convergence  of  x and  v occurs.  In  the  present  program, 
only  one  iteration  in  this  procedure  is  used,  but  the  error  introduced  by  this  simplification  is 
very  small.  Once  the  values  of  x are  known  from  the  beginning  of  the  path  to  the  various  ele- 
mental boundaries,  the  various  transmittances  are  simply  calculated  from  t(v)  = 1 - a(v)  using 
T . (2).  From  the  temperatures  of  the  various  path  elements,  I*c>nbe  evaluated,  and  the  inte- 
gration of  Eq.  (6)  performed. 

(U)  This  program  requires  a large  amount  of  storage,  and  its  execution  time  is  long,  so 
that  efforts  have  been  directed  toward  breaking  the  plume  into  elements  in  as  efficient  a manner 
as  possible.  In  some  regions  of  the  plume,  this  effort  has  led  to  a coarser  -.crementatlon  of  the 
plume  than  is  desirable.  From  experience,  it  has  been  determined  that  in  such  areas  the  contri- 
butions to  the  spectral-radiant  Intensity  are  somewhat  underestimated.  Further  work  is  needed 
to  completely  eliminate  this  bias. 

(U)  Further  details  of  the  operation  of  this  program  are  given  in  Appendix  IV. 


Appendix  IV 

DETAILS  OF  THE  PLUME- RADIATION  COMPUTER  CODE 

(U)  This  program  computes  the  radiant  intensity  of  the  plume  from  a rocket  or  jet  engine 
exhaust  as  observed  through  an  arbitrary  atmospheric  path.  The  elements  of  the  atmospheric 
path  are  computed  first  and  then  appended  in  front  of  the  elements  through  each  line  of  sight 
through  the  plume.  The  band-model  calculation  for  every  line  of  sight  is  thus  a nonisothermal 
band-model  calculation  including  the  path  and  that  line  of  sight. 

(U)  The  primary  coordinate  system  used  for  the  plume  calculation  is  a cylindrical  coordinate 
system  with  the  origin  at  the  center  of  the  nozzle  exit.  This  coordinate  system  is  shown  in  Fig.  6 
along  with  a generalized  plume  boundary  surface.  The  program  is  written  such  that  the  plume 
boundary  must  be  a conical  section  immediately  downstream  of  the  nozzle  exit,  followed  by  an 
optional  cylindrical  section,  This  surface  serves  only  to  define  the  region  over  which  integra- 
tion takes  place;  the  plume  parameters  can  be  arbitrarily  set  within  the  volume  enclosed  by  this 
surface.  This  primary  coordinate  system  is  denoted  by  subscript  1 in  the  figure.  The  axis  of 
this  system,  Zj,  is  coincident  with  the  plume  axis,  and  the  ^ = 0 direction  is  the  projection  in 
the  z j =0  plane  of  the  line  of  sight  directed  from  the  center  of  the  nozzle  exit  to  the  observer. 

(U)  The  second  coordinate  system  used  in  the  calculation  is  shown  in  Fig.  6 with  subscript 
2.  The  Zj,  axis  is  directed  along  the  line  of  sight  from  the  observer  to  the  plume.  Its  origin 
coincides  with  that  of  the  primary  coordinate  system.  The  <P2  = 0 direction  is  the  projection  in 
the  z2  = 0 plane  )f  the  plume  axis  directed  downstream  from  the  nozzle  exit.  The  integration 
of  the  plume  radiance  over  the  whole  plume  is  carried  out  by  projecting  the  plume  boundary  sur- 
face into  a plane  perpendicular  to  the  z2  axis. 

(U)  The  calculation  proceeds  as  follows: 

(1)  The  band-model  parameters  are  read  In. 

(2)  The  specifications  of  the  plume  boundary  (cone  angle,  nozzle  exit,  radius,  length  of 
conical  section  and  length  of  cylindrical  section)  are  read  In. 


(3)  The  observing  aspect  (in  degrees  from  nose-on  direction)  and  the  parameters  describ- 
ing the  atmospheric  path  between  the  plume  and  the  observer  (path  length,  altitudes  of 
both  the  observer  and  the  plume)  are  read  in.  Also,  the  required  atmospheric  proper- 
ties are  read  in  at  this  point. 

(4)  The  frequencies  for  which  calculations  are  desired  are  either  read  in  or  computed 
from  an  initial  value  and  increment. 

(5)  The  properties  of  the  plume  flow  field  (temperatures,  and  the  partial  pressures  of  the 
three  gases,  HjO,  C02  and  N2)  are  read  in  (in  spherical  coordinates  corresponding  to 
the  primary  coordinate  system).  An  option  is  also  available  in  which,  from  ambient 
and  core  values  read  in,  the  plume  properties  are  computed  according  to  an  internally 
specified  functional  relationship  within  the  plume  boundary.  If  this  option  is  used,  only 
core  and  ambient  values  for  temperature  and  the  three  partial  pressures  are  required. 

(6)  Various  control  parameters  are  read  in  which  specify  various  functions,  such  as  the 
fineness  of  the  integration  over  the  plume  (NMAX),  the  number  of  elements  into  which 
the  atmospheric  path  is  divided  (NA),  the  type  of  plume-data  read  in  (SETUPS),  and 
others  controlling  printing  and  plotting. 

(7)  The  temperature  range  over  which  the  plume  data  extend  is  read  in  and  then  divided 
into  a number  of  segments  (equal  to  NMAX)  according  to  the  2nd  power  of  the  tempera- 
ture. The  higher  the  temperature,  the  smaller  the  temperature  segment.  The  bound- 
aries of  the  segments  are  used  to  divtde  the  individual  lines  of  sight  into  elements. 

(8)  Subroutine  ATMO  is  used  to  divide  the  atmospheric  path  into  NA  elements.  At  this 
point,  the  length,  temperature  and  partial  pressures  of  HgO,  C02  and  N2  for  each  at- 
mospheric path  element  are  stored. 

(8)  Integration  over  the  plume  proceeds  as  follows.  A polar  element  in  the  plane  perpen- 
dicular to  the  z2  axis  is  constructed.  The  Initial  size  of  this  element  is  determined 
by  one  of  the  central  parameters,  NMAX.  A line  of  sight  through  the  center  of  this 
element  is  constructed  and  its  intersections  with  the  plume  boundary  surface  deter- 
mined in  the  z2,  r2,  and  <t>2  coordinate  system.  If  there  is  only  one  or  zero  points  of 
intersection,  the  current  line  of  sight  is  at  or  outside  the  plume  boundary.  The  pro- 
gram then  determines  where  this  line  of  sight  lies  with  respect  to  the  plume  boundary 
and  either  chooses  a new  element  or  terminates  the  computation,  whichever  Is  required. 

If  two  intersections  are  formed,  the  program  proceeds  to  divide  the  path  between  these 
two  Intersections  into  elements.  The  temperature  profile  along  this  line  of  sight  is 
examined  and  an  element  defined  whenever  the  temperature  profile  crosses  one  of  the 
boundaries  of  the  temperature  segments  defined  earlier.  The  length  of  each  of  these 
elements  is  stored  along  with  its  corresponding  mean  temperature  and  partial  pressures. 
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element  for  the  line  of  sight  currently  being  considered.  The  band-model  transmittance 
routine  TAU5  is  now  used  to  compute  elemental  transmittances  for  the  combination  of 
the  atmospheric  path  and  the  current  line  of  sight.  Radiance  values  at  the  observer 
are  now  calculated  and,  when  multiplied  by  the  area  of  the  current  element  projected 
into  a plane  perpendicular  to  t.ie  z2  axis,  are  stored  in  a running  sum  which  will  be 
the  plume  radiant  intensity  when  the  calculation  is  complete.  This  paragraph  is  re- 
peated for  all  frequencies  desired. 


(11)  An  average  temperature  for  the  current  line  of  sight  is  computed.  This  is  used  to 
modify  the  raHi-’l  size  of  the  next  element  in  the  integration.  The  next  element  is  now 
constructed  in  the  same  direction  but  incremented  radially.  The  radial  size  of  this 
element  is  based  on  the  average  temperature  of  the  previous  element  compared  to  the 
average  temperature  of  the  first  element  in  the  current  <t>2  direction.  The  first  two 
elements  in  any  <t>2  direction  have  the  same  radial  dimension,  but  the  smaller  the  aver- 
age temperature  of  succeeding  elements,  the  larger  the  radial  size  of  the  increment 
and  vice  versa.  The  value  of  <*>2  remains  constant  and  r2  is  incremented  until  a line  of 
sight  passes  outside  the  plume.  At  this  point  <f>2  is  incremented  and  elements  are  again 
constructed  starting  at  r2  = 0.  When  a new  element  has  been  defined,  the  calculation 
returns  to  the  point  of  paragraph  9 and  repeats  the  line-of-sight  and  radiance  calcula- 
tions. When  <t>2  has  been  incremented  over  its  full  range  (90°  in  the  case  of  90°  aspect, 
180°  in  the  case  of  all  others)  the  calculation  is  complete. 

(U)  The  running  sum  of  radiance  times  elemental  area  for  each  frequency  is  multiplied  by 
2 to  account  for  the  fact  that  integration  took  place  over  only  one  half  of  the  plume,  and  the  re- 
sults are  printed  as  the  radiant  intensity. 
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LIST  OF  SYMBOLS 


Constants  used  in  analytic  form  for  q 
average  line  spacing 


function  providing  interpolation  between  nearly 
weak  and  nearly  strong  approximations 


Ladenburg-Reiche  function 
spectrometer  slit  function 
absorption  coefficient 


averaged  absorption  coefficient  for  band  model 


2 -1 

spectral  radiance  W/ (cm  - sr-cm  ) 


2 -1 

black  body  spectral  radiance  W/ (cm  -sr-cm  ) 


,th 


black  body  spectral  radiance  of  i path  element 
spectral  radiance  computed  using  a band  model 
number  of  path  elements  in  sum 
line  strength 


equivalent  line  strength 

-1 

equivalent  width,  cm 


optical  depth  coordinate,  gm  cm 


total  optical  depth  of  path  of  interest,  gm  cm 


-2 


optical  depth  coordinate  of  it*'  path  element,  gm  cm 


dimensionless  optical  depth  for  Curtis-Godson 
approximation 


dimensionless  optical  depth  for  nearly  weak  - 
nearly  strong  approximation 


, 2tty 

line  overlap  parameter,  — -p- 


lina  overlap  parameter  averaged  over  the  path 
extending  from  the  observer  to  optical  depth  X 
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line  overlap  parameter  averaged  over  the  whole 
optical  depth  of  interest 
-1 

line  width,  cm 

line  width  averaged  over  the  path  extending  from 
the  observer  to  optical  depth  X 

line  width  averaged  over  the  whole  optical  "depth 
of  interest 

» 

interpolation  parameter  between  nearly  weak  - 
nearly  strong  approximation 
-1 

wavenumber,  cm 

wavenumber  of  line  center,  cm  1 
-1 

v-vq , cm 
transmittance 

averaged  transmittance  computed  using  a band  model 
averaged  transmittance  of  equivalent  non-overlapping 
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Introduction 

A recent  note  ^ was  concerned  with  the  application  of  a molecular 

(2  3) 

band  model  for  inhomogeneous  radiating  gases  * to  cases  in  which  a 
large  variation  in  temperature  is  encountered  along  the  optical  path.  The 
specific  problem  considered  was  the  calculation  of  thermal  radiation  from 
a relatively  small  combustion-product  source  viewed  through  a long 
atmospheric  path.  The  particular  band-model  used  in  that  study,  one 
previously  developed  at  this  laboratory  ’ , was  originally  formulated 

in  terras  of  two  options,  a "nearly-weak"  and  a "nearlv-strong"  line 
approximation,  with  a simple  procedure  for  interpolation  between  the 
two  for  intermediate  optical  depths.  Application  of  this  band  model 
in  its  original  form  yielded  anomalous  results;  indeed,  at  some  fre- 
quencies, the  calculated  values  of  apparent  radiance  with  atmospheric 
absorption  would  exceed  those  for  the  hot  source  alone.  The  reasons 
for  this  misbehavior  of  the  band  model  were  determined,  and  means  for 
Its  suppression  were  developed.  In  the  aforementioned  note,  some  pre- 
liminary results  of  the  use  of  the  modified  band  model  were  presented, 

i 

and  the  requirements  for  a more  complex  interpolation  procedure  was 
stated.  The  purpose  of  this  communication  is  to  present  the  details 
of  the  revised  formulation  of  the  model,  which  in  essence  represents 
a closer  approximation  to  the  exact  equation  of  transfer  than  the 
original  model,  or  comparable  ones  based  on  the  "Curtis-Godson" 
approximation. 
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The  rationale  in  the  development  to  follow  can  be  illustrated  by 
considering  the  growth  of  a single,  collision  broadened  spectral  line 
along  a non- iso thermal  path.  Once  the  formulation  for  such  a single 
line  is  developed,  it  is  carried  over  essentially  unchanged  to  the  case 
of  a random  band  of  overlapping  lines. 

The  solution  of  the  equation  of  transfer  describing  the  radiance  of 
a single  isolated  spectral  line  from  a general  non-iso thermal  source 
under  the  conditions  of  local  thermodynamic  equilibrium,  negligible 
scattering,  and  no  significant  source  of  radiation  behind  the  gas,  ie 


given  by; 


rxT  / 

x \ 

Vv) ' 

L j 

L*v(v,X)  k(v-vo,  X)  expl- 

o \ , 

k(v-vo,  X')dX'j  dX 

(1) 

X is  the  optical  depth  (mass  per  unit  cross  sectional  area)  coordinate  along 
the  line  of  sight  through  the  gas  and  L*v(v,X)  is  the  blackbody  spectral 
radiance.  The  latter  is  a function  of  frequency,  or  equivalently  wave- 
number,  v,  and  temperature,  T ■ T(X),  in  turn  a function  of  the  position 

along  the  line  of  sight  through  the  gas.  The  spectral  absorption  coefficient 

« 

of  the  single  line  being  considered,  k(v-VQ,  X)  ia  a function  of  the  spectral 
distance  from  the  line  center  v-vq,  and  the  optical  depth  coordinate  along 
the  path  X.  Ly(v)  is  the  spectral  radiance  of  the  gaseous  source  at  wave- 
number  v due  to  the  single  spectral  line.  X^  is  the  total  optical  depth  — 

of  the  gaseous  source.  The  total  radiance  due  to  the  line  is  obtained  by 
integrating  over  all  frequencies,  thus 


r*L 


fX 


L*v(v,X)  k (v-vq,  X)  exp [ - 


k(v-vo,  X')dX'|  dX  d(v-vft)  (2) 
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By  introducing  the  transmittance  defined  as, 

r 

rX 


t(v-v  , X)  = exp 


k(v-v  X')  dX' 


(3) 


interchanging  the  order  of  integration  and  assuming  that  L*  (v,X)  does  not 
vary  significantly  over  the  frequency  region  involved,  the  equation  for  L 
can  be  written: 


ft. 

»oo 

L ■ - 

LW  « dX 

t(v-vq,  X)d(v-v‘o) 

. 

o 

J 

dX 


(Z  a) 


The  equivalent  width,  W,  is  defined  as 
W(X)  « 1 - 

so  that  the  expression  for  L becomes 


t(v-v  , X)  d(v-v  ) 
o’  o 


(V) 


L - 


dW 

L\<V  X>  f dX 


(£4) 


Since  W Is  only  a function  of  X,  the  above  equation  can  be  written 

. _ f W(V 


L*v(vo,  X)  dW(X) 


(5) 


By  dividing  up  the  path  into  N approximately  isothermal  elements,  the 
expression  for  L can  be  written: 


N 

L - E L*  [W(X  ) - W (X.  )]  v 

i»l  V»1  1 1-1 


(6) 


where  X.^  is  the  optical  path  from  the  observer  to  the  far  boundary  of  the 
i isothermal  element.  W(X^)  is  the  equivalent  width  corresponding  to  that 
path. 
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calculatln8  «‘e  vlue.  of  WfX^  it  1,  that,  in  generel.for  eon- 

lsochemel  paths,  the  lntegtala  over  X and  v-vo  are  not  separable.  To  allow  ' 
such  separation,  en  epproximetion  la  usually  made  concerning  the  variation 
of  the  line  .hep.  along  the  peth.  The  Curti.-Godson  (4-5'6>  approaleatlon 
i.  .oet  commonly  used  for  thi.  purpoae.  toother,  that  i,  „.ed  by  these 
author.  (ll3),  1.  the  "nearly-we.k,  ne.rly-strong"  line  approximation,  here- 
after referred  to  a.  the  Ntf-HS  approximation.  Both  approximation,  have  as 
their  purpose  the  seperatlon  of  the  integral,  ovar  X and  v-vQ  and  the  re- 
placement of  the  non-iso  thermal  path  with  an  equivalent  isothermel  path. 

Both  yield  **tla factory  value,  of  V(X)  for  many  cases,  especially  where  the 
temperature  and  concentration  variations  are  not  too  large. 

However,  in  the  evaluation  of  eq.(6>,  the  accuracy  of  the  integration 
depends  not  so  much  on  the  accuracy  of  the  values  of  W(X),  but  on  the  accuracy 
of  differences  between  values  of  W(X),  i.e.,  on  the  derivative  of  W with 
respect  to  X.  Differentiations  of  approximations  are  risky  and  indeed  this 
is  the  point  at  which  the  anomaly  mentioned  earlier  appears.  An  expression 
which  is  a reasonable  approximation  for  W will  not  always  yield  a reasonable 
approximation  to  dW/dX.  To  illustrate  this,  let  us  compute  the  derivative  of 
W with  respect  to  X using  each  of  the  two  approximations  mentioned  earlier; 
the  Curtis-Godson  approximation,  and  the  NV-NS  approximation. 

Before  this  can  be  accomplished  it  will  be  necessary  to  obtain  an  expression 
for  the  equivalent  width  of  a collision  broadened  line  in  an  isothermal  path. 

The  absorption  coefficient  in  this  case  is: 


k(v-vo,  X') 


S(X*) 

x 


-T(X-) 


yZ(X')  + <v-v  )2 


(7) 
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where  S is  the  line  strength  and  y its  half  width.  Upon  substitution  into 

/ ' 

eq,  (3)  specialized  to  an  isothermal  path,  and  with  the  use  of  eq. 

after  some  reduction,  the  following  expression  for  W is  obtained  ^ 

W(X)  - 2 try  f(x)  (7a) 

where  f is  the  Ladenburg-Reiche  function 

f(x)  = x e x [J  (i  x)  -*i  J.(i  x)]  C7b) 

a 1 

sx 

aad  x " 2tty  is  the  dimensionle88  optical  depth.  Jq  and  J.^  are  Bessel 
functions  of  the  first  kind,  of  order  0 and  1 respectively.  Tables  of  f(x) 
are  found  in  the  literature 

Curtis-Godson  Approximation 

The  Curtis-Godson  approximation  has  been  defined  in  several  ways.  The 
form  described  by  GOODY ^ is  used  here.  It  is  a method  whereby  the  equivalent 
width  of  an  absorption  line  as  viewed  through  a non -iso thermal  path  can  be 
approximated  by  the  equivalent  width  through  some  corresponding  isothermal 
path.  The  approximation  is  made  by  defining  an  equivalent  average  half  width 
and  line  strength  for  the  path  extending  from  the  observer  to  some  point  X in 

i 

the  path. 


and 


(8) 


This  is  tantamount  to  saying  that  the  denominator  in  the  non-isothermal 
absorption  coefficient,  eq.  (7),  can  be  taken  to  be  independent  of  X^  and 


B 


to  have  the  form,  y 2 (x)  + (v-v  )2,  so  that 

e O 

[X  , . „ Y.(X)  S (X)  X 

k(v-v  , x ) dX"  - ~ £■ 

J ^ O _ 9 


TIlYe2(X)  + (v-vQ)2] 


From  this,  it  follows  that 


W(X)  - 27Tye(X)  f(x) 

where  x,  the  dimensionless  optical  depth,  is  now  defined 


x(X)  = 


S(x')dX' 

2irve(X) 


Eq.  (10)  can  now  be  used  to  investigate  the  behavior  of  dW/dX  in  a non- 
i6o thermal  gas.  Differentiating  (10)  and  making  use  of  (8)  yields  after 
some  rearrangement: 


I ■ s<« 


i *■  \ ytm) + » - M} 


By  Mapping  th.  variation  with  optical  depth  X Into  a variation  with  at 

less  optical  depth,  x,  this  can  be  rewritten 

(*-£«)  (13) 

We  can  now  place  physical  limitations  on  [1/S(X)] (dW/dX) . First  of  all,  at 
X - 0,  the  quantity  (1/S (X)] (dW/dX)  must  be  unity  corresponding  to  the  region 
of  linear  growth  for  small  optical  depths.  Secondly,  it  must  be  less  than 
one  at  larger  optical  depths  as  the  square  root  region  is  approached.  Thirdly, 
it  must  remain  non-negative  for  all  x so  that  the  equivalent  width  continues 


f ' 1 


<J>  jfl 


to  grow  with  increasing  x.  The  variation  of  [1/S(X) ] (dW/dX)  as  given  by 

eq.  (7)  is  shown  in  Figure  1 as  a function  of  x with  y/yb  as  a parameter. 

Note  that  Y/Ye  is  the  ratio  of  the  half  width  at  x to  its  average  value 

taken  over  all  the  path  before  it.  For  y/y  < 2 the  function  is  well 

e 

behaved,  but  it  exhibits  anamalous  behavior  for  y/y^  >2.  In  the  latter 
case  the  function  becomes  larger  than  1 for  intermediate  values  of  x, 
thus  indicating  a rate  of  increase  of  equivalent  width  greater  than  that 
in  the  linear  region.  This  behavior  is  a manifestation  of  the  errors 
involved  in  the  Curtis-Godson  approximation. 

These  errors  are  mentioned  by  GOODY  ^ ; and  are  treated  in  detail  by 
DRAYSON^5\  However  their  analyses  dealt  only  with  the  absorption  (or 
equivalent  width).  More  related  to  this  study  is  the  work  of  WALSHAW  and 
RODGERS who  analyzed  the  effect  of  the  Curtis-Godson  approximation  on 
the  derivative  of  transmittance  with  optical  depth,  for  several  band  models. 

Let  us  now  examine  the  behavior  of  the  NW-NS  approximation  in  a similar 
fashion. 


Nearly-Weak  Nearly-Strone  Line  Approximation 


This  approximation  is  similar  to  the  Curtis-Godson  approximation  in 

that  an  average  for  the  line  half  width  is  defined.  However  this  average 

is  fixed  for  the  path  in  question  and  is  defined  as  the  average  over  the 

total  path  X^ . Thus  ;XT 

L L Ye(X')  S(X')dX' 

Y - ' (14) 

eL 


(*1 

S(X') 

r\ 


Then  the  dimensionless  optical  depth  is  defined  as: 


3 


Efc, 


x(X)  - 


2iry 


eL 


S(X')  F(X')dX' 


(15) 


F(X')  is  a function  which  must  be  unity  in  the  nearly-weak  approximation 

(ID 


and  y/y  In  the  nearly-strong  approximation 
6 L 


A form  must  be  chosen 

for  the  variation  of  F(X')  so  that  it  produces  a reasonable  transition  bet-* 

» 

ween  the  two  approximations  for  intermediate  optical  depths.  The  equivalent 
width  is  given  again  as 


W(X)  - 2itYeLf(x) 


(16) 


Taking  the  derivative  of  eq.  (16)  we  obtain,  for  the  nearly  weak  approximation 

(17) 


dW 

dX 


S(X)  - 

ax 


and  for  the  nearly  strong  approximation 


dW 


dX 


- S(X) 


'eL 


dx 


(18) 


[1/S(X)  ] (dW/dX)  according  to  eqs.  (17)  and  (18)  is  shown  in  Figure  2.  It 
can  be  seen  that  the  nearly-weak  approximation  never  becomes  physically 
unrealistic,  i.e.,  [ 1/S (X)] (dW/dX)  <_  1;  however,  for  the  nearly-strong 
approximation,  the  quantity  [1/S (X)] (dW/dX)  becomes  greater  than  one  at 


smaller  values  of  x for  > 1.  This  difficulty  is  accounted  for  in 


'eL 


practice  by  the  proper  choice  of  F(X').  However,  we  are  left  with  no 
information  upon  which  to  base  F(X')  at  this  stage. 
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It  appears  that  neither  the  NW-NS  approximation  (unless  F(X)  is  correctly 
and  accurately  specified)  nor  the  Curtis-Godson  approximation  give  suffi- 
ciently accurate  values  for  dW/dX  over  the  full  range  of  optical  depths  and 
for  all  possible  values  of  y/y&»  The  treatment  of  the  previous  paragraphs 
serves  to  point  out  those  regions  where  the  approximations  are  not  only 
inaccurate  but  yield  physically  meaningless  results  as  well.  On  the  other 
hand,  it  should  be  noted  that  the  anomalies  and  errors  in  the  curves  for 
Cl/sDfclW/d^as  given  by  the  Curtis-Godson  approximation  are  not  serious  as 
long  as  Y /Y  never  exceeds  a value  of  about  2.  This  condition  has  been 
shown  to  be  satisfied  ^ *5»10)  for  most  atmospheric  paths  and  only  becomes 
seriously  violated  for  highly  non-homo geneous  paths.  However,  as  an  ex- 
ample  of  the  latter  taken  from  our  studies  ’ , hot  products  of  combustion 

viewed  through  long  atmospheric  paths  yield  values  of  y/y^  as  large  as  20. 

To  obtain  a ^more  acceptable  expression  for  dW/dX,  consider  the  growth 
of  an  isolated  spectral  line.  The  equivalent  width  is  given  exactly  by 


W(X) 


■ L exp  [•!< 


k (Av,  X')  dXjd(Av) 


1 


(19) 


and  its  derivative,  in  terns  of  a Lorentz  profile,  by 


dW 

dX 


i r scx)v(x) 

* y2(X)  + (Av)2 


SpxQQ.  dX1*]  d (Av) 
Y (X'>  + (Av)z 

m 


(20) 


where  Av  - v - vq.  In  the  original  treatment  of  the  Curtis-Godson  and  the 
NW-NS  approximations,  the  separation  of  frequency-dependent  and  path- 


dependent  variables  was  accomplished  by  substitution  of  an  effective  value, 
Yft>  in  both  parts  of  the  integrand.  At  this  point  we  will  use  such  an 
approximation  only  in  the  exponential  term.  This  is  the  key  step  in  this 
derivation  and  is  the  only  point  of  difference  from  previous  treatments. 


Thus,  in  teras  of  Yfi  - Yg(X)  defined  by  equation  8 and  x by  eq.  (11), 
eq.  (20)  becomes 


dW  _ 1 f 

dX  " it  ] 


.S(x)y.&L Up  -2x  Ye00 

Y (X)  + (Av)2  |_ 


YeCX) 

Y#2(X)  (Av)2 


d(Av) 


Note  that  this  is  equivalent  to  making  a Curtis-Godson  substitution  in  the 

expression  for  the  derivative  of  W rather  than  in  the  expression  for  U 

» 

itself.  Thus  the  process  of  taking  a derivative  after  making  an  approxi- 
mation has  been  replaced  by  one  in  which  the  derivative  is  taken  first  and 
then  the  approximation  made.  Hence  this  treatment  should  be  inherently 
better  for  calculating  radiances  than  the  techniques  described  earlies. 
Another  way  of  comparing  this  approximation  to  the  previous  two  is  by 
considering  dW/dX  aa  the  limiting  value  of  the  contribution  to  W from  a 
small  element  of  path,  Ax,  viewed  through  the  path  ahead  of  it,  X.  In 


(21) 


the  previous  two  approximations  the  line  shape  of  the  element  AX  and  the 
path  in  front  of  It  are  both  given  by  the  line  shape  for  the  equivalent 
homogeneous  path.  In  the  current  approximation  the  element  AX  is  taken  to 
have  its  true  line  shape  while  the  path  X in  front  of  it  is  taken  to  have 
an  equivalent  homogeneous  shape.  Eq.  (21)  has  been  checked  against  the 
exact  relation,  eq.  (20),  for  a two  layer  non-isothermal  path  and  found 
to  be  very  accurate.  This  is  reasonable  based  on  the  physical  interpreta- 
tion given  above.  Further  investigation  of  its  accuracy  is  to  be  the  sub- 
ject of  future  work.  It  would  be  possible  t<J  go  on  from  here  and  fully 
develop  a procedure  for  treating  isolated  Lorentz  lines,  but  since  our 
interest  is  ultimately  in  a band  model,  which  treats  the  average  of  many 


lines,  let  us  turn  our  attention  at  this  point  to  a random  band  model, 
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treating  it  in  a manner  analogous  to  that  used  above  for  the  isolated  line. 
Band  Model  Considerations 

The  calculation  of  an  average  spectral  radiance  at  the  boundary  of  a 
radiating  body  of  gas,  in  a frequency  interval  encompassing  a number  of 
rotational  lines,  is  represented  by  the  relation 

V f s(v'-v)  Lv(v')  dv'  (22) 
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where  g(v,v")  is  an  averaging  function,  e.g.  the  slit  function  of  a spectro- 
meter of  moderate  resolution,  such  as  those  which  produce  original  absorption 
spectra  from  which  band-model  parameters  are  extracted.  The  exact  spectral 
radiance  Ly  is  given  by  the  equation  of  transfer,  eq.  (1)  and  (3),  so  that 
eq.  (22)  can  be  written 
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where  t iB  the  perfectly  resolved  transmittance,  and  the  primes  denote 
variables  of  integration.  • 

The  order  of  integration  in  eq.  (23)  can  be  inverted,  and  for  frequency 
intervals  small  enough  that  the  Planck  function  L^*(T)  is  essentially  constant, 
the  result  is  the  band-model  expression: 
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where 
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t(v)  = 

Thus  7 is  seen  to  be  the  "spectral"  transmittance  measured  by  a conventional 
instrument  or  calculated  by  a band-model. 

The  evaluation  of  eq.  (25)  to  obtain  the  average  spectral  radiance 
is  done  by  first  computing  t at  a number  of  pointa  along  the  path.  The 
blackbody  radiance  L^*  is  then  represented  as  a function  of  t,  and  Integrated 
numerically  over  t,  using , for  instance,  a simple  summation.  The  various 
band-models  in  current  use,  including  those  based  on  the  Curtis-Godson  approx- 
imation and  the  preaent  model  in  lta  original  form,  do  yield  values  for  the 
transmittance  with  an  accuracy  sufficient  for  most  engineering  applications. 
However,  as  in  the  case  of  isolated  lines,  the  use  of  eq.  (25)  in  essence 
involves  a differentiation  of  the  expression  for  7 with  respect  to  the 
optical  path.  Differentiation  of  course  accentuates  errors  and  uncertainties; 
in  the  isolated  line  case  presented  above,  it  was  the  derivative  of  the 
equivalent  width  that  exhibited  the  misbehavior.  Therefore,  analogously  to 
the  case  for  isolated  lines,  an  expression  which  yields  reasonable  values 
forTwill  not  necessarily  yield  reasonable  values  for  dT  and  hence  Ly. 

What  is  needed  is  an  expression  for  7 having  a derivative  which  will  yield 
more  realistic  values  of  the  difference  in  transmittance  with  a specified 
incremental  increase  in  physical  path.  Such^  an  expression  can  be  developed 
as  follows. 

The  transmittance  for  a band  of  spectral  lines  with  some  overlspping  can 
(4) 
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be  expressed  as 


where  Tn  is  the  transmittance  for  an  equivalent  band  of  well-isolated  non 
overlapping  lines.  Differentiating  eq.  (26)  with  respect  to  the  optical 


depth  X yields 

I 

di_  — di 

dX  "T  (27) 

dX 

which  indicates  that,  with  a reasonably  accurate  value  for  7,  the  accuracy 
of  dx/dX  depends  directly  on  that  of  d7n/dX.  Accordingly,  an  improved 
relation  for  the  latter  will  be  sought. 


In  our  current  work  we  have  used  the  NW-NS  approximation  exclusively  and 
our  consideration  will  be  limited  at  this  point  to  this  model,  seeking  a 
specification  of  the  function  F(x')  which  will  yield  accurate  values  of  a 
transmittance  derivative.  The  Curtis-Godson  approximation  contains  no 
unspecified  parameter  which  can  be  so  chosen,  so  it  will  not  be  considered 
further.  For  the  NW-NS  approximation  (11) 


Tn  " eeL  f(x) 


where  the  argument  of  the  Ladenburg  Reiche  function  is  in  this  case 
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n(x') 


in  which  k - k(v,T)  is  the  first  band-model  parameter,  the  average  absorption 
coefficient,  identifiable  as  the  average  line  strength  to  spacing  ratio,  and 
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0 •«  0 (v,T)  is  the  second  band-model  parameter,  the  line  overlap  factor, 
identifiable  as  2ir  times  the* average  line  width  to  spacing  ratio  . The 
effective  value  0 is  defined  analogously  to  for  the  single  case: 


non  , 

F(X')  has  here  been  replaced  by  II ^X-  , a suitable  functional  form 

eL 

for  F(X').  Here  n is  an  interpolation  parameter  between  the  nearly  weak 
and  nearly  strong  approximations,  varing  from  0 for  the  nearly  weak  to  1 
for  the  nearly  strong. 

In  an  earlier  study  , an  investigation  of  the  form  for  n is 
described.  As  a result  of  that  study  several  empirical  forms  for  n 
were  determined  which  appeared  to  be  reasonably  uaeful.  However,  the 
information  used  to  make  these  determinations  was  not  greatly  sensitive 
to  values  of  n and  hence  accurate  information  about  its  form  could  not 
be  extracted.  One  of  these  empirical  forms  is: 
n - 0 for  x < 1 
n - (x  - l)/9  for  l x <_  10 
n ■ 1 for  x > 10 

This  appears  to  make  n a function  of  x rather  than  X'*,  and  the  determination 
of  x an  interative  procedure.  However,  in  practice  the  integral  for  x is 
evaluated  by  a summation  and  a value  for  n for  a particular  term  is  deter- 
mined from  the  value  for  x calculated  from  the  sum  of  the  preceding  terms. 
Thus  n is  uniquely  defined  for  each  X*. 

Differentiating  eq.  (28)  end  utilising  eq.  (30)  yields 
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(31) 
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Which  shows  that  the  derivative  of  dT^dX  varies  with  the  absorption 

coefficient,  as  would  be  expected;  the  derivative  of  the  function  f(x), 

and  the  local  values  of  0(X)  and  n(X).  The  physical  significance  of  eq. 

(31)  can  be  seen  in  Figure  3,  a normalized  plot  of  -(1/k)  d7  /dX  vs  x 

n 

with  0/egL  as  a parameter.  In  the  limit,  as  x+0,  the  band  model  reduces, 

as  it  should,  to  a linear  growth  law,  so  that -(1/k)  dT/dX  approaches  unity. 

I 

However,  this  function  must  be  less  than  unity  for  all  values  of  X greater 
than  zero,  in  accordance  with  the  physical  requirement  that  the  trans- 
mittance decrease  monotonically  with  optical  depth.  In  Figure  3,  the 
empirical  specifications  for  n lead  to  a violation  of  this  constraint  for 
certain  values  of  3/ and  x.  This  misbehavior  is  analogous  to  that 
shown  earlier  for  isolated  line  growth.  We  now  search  for  a better 
approximation  to  dt^/dx.  This  quantity  is  given  by 

dx  — 

n 1 dW 

dx  d dx  • (^2) 

where  W is  an  average  equivalent  width  and  d is  an  average  line  spacing  in 
the  spectral  interval  under  consideration.  For  lines  of  equal  strength 
and  width,  after  some  rearrangement,  eqs.  Q2)  and  (21)  yield: 
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where  x denotes  the  dimensionless  optical  depth  defined  analogously  to  x as 
defined  In  eq.  (11)  and  used  In  eq.  (21).  Rewritten  In  terms  of  the  two 


bandmodel  parameters,  It  becomes 

f k (X')  dX" 

x ■ ^o  (34) 

eeoo 


where  @£(X)  Is  the  average  value  of  the  second' band  model  parameter  over  the 

path  from  the  observer  to  the  point  X,  analogous  to  the  definition  of  v (X)  In 

e 

•q.  (8) 
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The  distinction  between  0eL  (a  constant,  since  it  is  determined  over  the  whole 
path),  defined  in  eq.  (30)  and  0£(X)  (a  function  of  X,  the  optical  depth  coor- 
dinate) should  be  noted. 

The  change  of  variablee  tan  6/2  ■ 2n  Av/0ed  and  a trlgnc, metric  substitution 
yields  after  some  manipulation 
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Eq.  (36)  haa  been  integrated  numerically;  the  results  are  plotted  in  Figure  4, 
again  as  -(1/k)  dTn/dX  vs  x with  0(X)/0e(X)  a?  a parameter.  The  resulting  curves 
are  well  behaved  and  exhibit  none  of  the  anomalies  seen  in  Figure  3.  Notice  that 
the  curve  for  0/0^  "lie  identical  in  both  figures. 

The  interpolation  parameter  n can  now  be  evaluated  by  equating  eqs.  (36)  and 
(31).  Thesa  two  equations  involve  different  sets  of  variables,  i.a.  eq.  (36)  is 
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The  relation- 


in  terms  of  x and  0 (X),  and  eq.  (31)  is  in  terms  of  x and  0 

6 cL 

ship  between  x and  x is  complex,  depending  on  q,0  (X)  and  0 in  an  integral 

e eL 

form.  Furthermore,  whereas  x is  a unique  function  of  the  optical  path  extend- 
ing from  0 to  X,  x is  not  a unique  function  in  the  same  sense,  depending, 
through  0^,  on  the  path  between  X and  X^  as  well.  Such  dependence  is  not 
physically  realistic  but  appears  as  an  integral  part  of  the  NW-NS  approxima- 

i 

tlon.  In  determining  a form  for  q this  non-uniqueness  must  be  adequately 
accounted  for  so  that  the  computed  value  of  (1/k)  (dx  /dX)  is  unique  despite 
the  non-unique  character  of  x.  It  is  very  difficult  to  insure  the  existence  of 
the  latter  condition  except  for  certain  types  of  paths.  In  particular  if  0 dif- 
fers significantly  from  its  average  value,  0 , only  in  a very  small  portion  of 

the  optical  path  then  the  contribution  from  that  portion  of  the  path  to  the  inte- 
gral defining  x is  negligible  compared  to  the  contributions  from  the  remainder  of 
the  path.  In  this  case  the  following  approximate  relation  holds  i 


The  above  restriction  also  implies  that 
0 (x)  a constant  - 0 


and  hence  that  v 

x a x 

Thus,  by  restricting  ouraslves  to  cases  where  0 (X)  is  not  a strong  function  of 
X and  allowing  large  variations  in  0(X)  only  in  very  small  fractions  of  the  opti- 
cal path,  values  for  n can  be  found  in  terms  of  x and  0/0  (except  for  0/0  - 1 

CL  6L  / 
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where  n is  indeterminate).  The  resultant  values  of  n so  obtained,  as  a function 

I ' 

of  x and  8/8  , are  shown  in  Figure  5;  using  these  values  in  eq.  (31)  will  in- 
eL 

sure  that  the  derivative  dr/dX  will  equal  that  given  by  eq.  (36).  It  should  be 
noted  that  x is  computed  by  a sum  so  that  the  value  of  n valid  for  a given  term 
of  the  sum  is  determined  from  the  sum  of  the  preceding  terms.  For  convenience 
in  application,  simple  analytical  approximations  were  made  to  the  curves  of  Fig. 
5 by  fitting  them  with  functions  of  the  form 

n(x,  8/8  ) " * f A 


(37) 


x + B (8/6,) 


n(x,  8/8  ) - X C (6/6e)- 


E (8/8  ) + at  D CB/80) 


(38) 


For  any  particular  Instance,  the  choice  between  eq.  (37)  and  (38)  is  made  on 
the  basis  of  the  quantity  s 


E(8/B,) 


1 + 0.185  6/8. 


(39) 


If  x > E,  eq.  (37)  is  used;  otherwise  eq.  (38)  is  used.  .Values  for  A,  B,  C, 

and  D are  tabulated  in  Table  I.  Values  of  Ac  /d X obtained  by  use  of  Table  I 

n 

in  place  of  the  values  shown  for  n in  Figure  5 are  sufficiently  accurate  for 
swat  applications  of  this  band-model. 

The  limitation  on  the  variation  of  6(X)  end  6e(X)  mentioned  earlier  is  not 
as  severe  as  it  might  first  appear,  being  satisfied  in  many  cases  of  interest. 
For  instance,  in  cases  where  a relatively  small  hot  gaseous  source  is  to  be 
viewed  through  a long  cool  path,  the  condition  is  approximately  met.  In  such 
a case,  the  value  of  6 is  almost  invariant  over  the  whole  atmospheric  path 
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which  contains  the  largest  part  of  the  optical  path.  0 varies  substantially 

/ 

only  in  the  relatively  snail  portion  of  the  path  comprising  the  hot  gaseous 
source.  Hence  the  condition  is  approximately  satisfied. 

It  should  be  noted  that  the  introduction  of  an  interpolation  parameter  is 
somewhat  of  an  ad  hoc  procedure  having  little  physical  basis  and  is  not  the 
only  way  eq.  (36)  could  be  introduced  into  a computational  procedure.  An  al- 
ternative way  to  use  the  above  results  would  be  the  direct  application  of  eq. 

(36)  into  a band-model  computer  code.  The  resultant  values  of  dT  /dX  could 

n 

then  be  used  in  eq.  (27),  with  in  eq.  (26)  now  evaluated  in  terms  of  x and 
6e(X).  Such  a procedure  would  obviate  the  need  for  a distinction  between  the 
nearly-weak  and  nearly-strong  approximations  and  hence  an  interpolation  proce- 
dure. Unfortunately,  in  the  use  of  eq.  (36),  an  analytic  form  for  the  solution 
of  the  integral  is  not  easily  obtained.  Although  a solution  has  been  obtained 
in  the  form  of  an  infinite  series,  this  represents  little  improvement  over  a 
table  of  numerical  values,  Whether  there  would  be  a net  gain  in  computational 
efficiency  remains  to  be  seen.  However,  it  might  be  expected  to  be  somewhat 
more  satisfactory  in  that  the  restriction  on  the  variation  of  6 (X)  compared  to 

t 

required  in  the  evaluation  of  q would  not  be  needed.  This  alternative  mode 
of  the  band-model  formulation  is  one  of  the  subjects  of  a continuing  investiga- 
tion. 

In  summary,  we  have  shown  how  both  the  Curtis-Godson  approximation  and  the 
nearly-weak, nearly-strong  approximation  can  introduce  anomalies  into  the  radi- 
ance computed  over  a highly  non-isothermal  path  for  an  isolated  spectral  line. 
We  have  than  applied  a Curtis  Godson  type  approximation  to  the  expression  for 
the  equivalent  width  gradient  along  the  path,  producing  a formulation  which 
yields  greatly  Improved  radiance  values,  exhibiting  no  anomalies.  Generalizing 


125 


to  a random  band  model,  the  expression  concerning  the  equivalent  width  gradient 

* / ' 

has  been  used  to  obtain  an  interpolation  between  the  nearly  weak  and  nearly 
strong  approximations  valid  for  that  band  model. 

The  most  important  result  presented  herein  is  the  form  of  the  transmittance 
gradient  given  in  eq.  (36).  It  is  the  basis  of  a computational  procedure  for 
computing  the  radiance  from  highly  non -iso thermal  gaseous  sources  which  yields 
significant  improvements  over  previously  used  ‘techniques.  It  has  been  shown  to 
be  a better  approximation  than  the  more  commonly  used  approximations,  and  yields 
physically  realistic  results  for  all  values  of  the  path  and  band  model  parameters. 
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